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Abstract 

This  research  extends  the  theory  and  understanding  of  the  laser  speckle  imag¬ 
ing  technique.  This  non-traditional  imaging  technique  may  be  employed  to  improve 
space  situational  awareness  and  image  deep  space  objects  from  a  ground-based  sensor 
system.  The  use  of  this  technique  is  motivated  by  the  ability  to  overcome  aperture 
size  limitations  and  the  distortion  effects  from  Earth’s  atmosphere. 

Laser  speckle  imaging  is  a  lcnsless,  coherent  method  for  forming  two-dimensional 
images  from  their  autocorrelation  functions.  Phase  retrieval  from  autocorrelation  data 
is  an  ill-posed  problem  where  multiple  solutions  exist.  This  research  introduces  po¬ 
larization  diversity  as  a  method  for  obtaining  additional  information  so  the  structure 
of  the  object  being  reconstructed  can  be  improved.  Results  show  that  in  some  cases 
the  images  restored  using  polarization  diversity  are  superior  to  those  reconstructed 
without  it. 

This  research  presents  statistical  analysis  of  the  observed  data,  two  distinct 
image  recovery  algorithms,  and  a  Cramer-Rao  Lower  Bound  on  resolution.  A  math¬ 
ematical  proof  is  provided  to  demonstrate  the  statistical  properties  of  the  observed, 
noisy  autocorrelation  data.  The  algorithms  are  constructed  using  the  Expectation- 
Maximization  approach  and  a  polarization  parameter  that  relates  two  independently 
observed  data  channels.  The  algorithms  are  validated  with  computer  simulation  and 
laboratory  experiment.  Comparison  is  made  to  an  existing  phase-retrieval  technique. 
The  theoretical  lower  bound  is  developed  for  comparing  theoretical  performance  with 
and  without  polarization  diversity.  The  results  demonstrate  the  laser  speckle  imaging 
technique  is  improved  with  polarization  diversity. 
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Statistical  Image  Recovery  From  Laser  Speckle  Patterns 


With  Polarization  Diversity 


I.  Introduction 

Synthetic  aperture  Laser  Detection  and  Ranging  (LADAR)  imaging  is  investi¬ 
gated  in  this  research.  This  research  provides  an  improved  phase  retrieval  algorithm 
for  recovering  remotely  sensed  objects  from  non-imaged,  laser-speckle,  intensity  data. 
The  innovative  aspects  investigated  as  part  of  this  effort  include  the  addition  of  po¬ 
larization  diversity  and  a  statistical  approach  to  object  recovery  that  improves  upon 
existing  techniques.  Improvement  is  demonstrated  via  computer  simulation,  labora¬ 
tory  experiment,  and  a  theoretical  bound  on  resolution.  This  work  is  a  continuation 
of  previous  research  efforts  [3,31-33]. 

Fundamentally,  the  Department  of  Defense  desires  cost  savings  and  improved 
image  resolution  beyond  what  can  be  achieved  with  existing  optical  reconnaissance 
systems.  In  general,  with  larger  apertures,  greater  diffraction-limited,  resolution  is 
achieved.  There  are  applications,  such  as  space-borne  sensors,  where  optical  imaging 
is  preferred;  however,  the  size  of  the  optic  is  limited  by  weight  and  cost  constraints. 
Practically,  for  a  satellite  system,  weight  and  cost  savings  can  be  achieved  at  the 
expense  of  power  and  computer  processing.  Weight  and  cost  savings  can  potentially 
be  achieved  with  large,  synthetic  apertures  without  a  focusing  lens.  Large,  synthetic 
apertures  designed  with  many  light-weight  and  relatively  inexpensive  detector  ele¬ 
ments  save  production  expense  and  weight  compared  with  large,  monolithic  optical 
elements  found  in  traditional  optical  systems.  Ideally,  an  even  larger  overall  aperture 
can  be  synthesized  with  many  small  elements  as  compared  to  a  cost  equivalent  optical 
system. 
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Along  with  cost  savings,  the  goal  is  to  improve  upon  available  image  resolution 
and  overcome  detection  and  atmospheric  noise.  For  ground-based  sensors  looking  sky¬ 
ward,  large  optical  systems  are  developed;  however,  the  Earth’s  atmosphere  severely 
limits  performance.  It  will  be  shown  in  this  dissertation  that  laser  speckle  imaging 
without  optics  in  ground  systems  avoids  the  effect  of  imaging  through  atmospheric 
turbulence  at  the  aperture.  Turbulence  effects  are  not  observed  because  the  per¬ 
turbed  phase  information  is  completely  destroyed  in  the  collection  process  [3,44], 
Only  scintillation  effects  due  to  layered  atmosphere  above  the  aperture  plane  need  be 
considered  [44].  However,  without  the  object  phase  information  we  must  estimate  the 
phase  leading  to  the  classical  phase  retrieval  problem  discussed  later  in  this  document. 

This  chapter  introduces  the  dissertation  research  and  its  documentation.  The 
operational  motivation  for  conducting  the  research  is  first  provided  in  Section  1.11 
Technical  motivation  for  this  research  is  provided  in  Section  1121  This  is  followed  by 
a  summary  of  Research  Contributions  in  Section  11.31  The  chapter  concludes  with  a 
Dissertation  Overview  in  Section  11741 

1.1  Operational  Motivation 

Accurate  and  reliable  Space  Situational  Awareness  (SSA)  is  an  important  re¬ 
quirement  for  the  United  States  Air  Force.  General  C.  Robert  Kehler,  Commander, 
Air  Force  Space  Command,  in  a  speech  to  the  Air  Force  Association  in  November 
2007  said, 


“...a  critical  important  thing  for  us  and  that  is  to  get  better  at  space 
situational  awareness.  That’s  one  of  our  top  priorities  in  the  command. 

It’s  going  to  remain  one  of  our  top  priorities  on  my  watch.  And  leaders 
flat  simply  have  got  to  get  better  about  knowing  what’s  up  there,  tracking 
what’s  up  there,  understanding  the  intent  of  things  that  are  up  there,  and 
knowing  those  pieces  in  real  time.”  [27] 

The  space  environment  continues  to  be  more  complex  due  to  increased  on-orbit  sys¬ 
tems  as  well  as  debris.  The  United  States  must  monitor  all  on-orbit  systems  and  debris 
for  safeguarding  and  positive  control.  The  SSA  mission  demands  timely  knowledge 
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of  debris  and  orbiting  systems  over  a  vast  region  of  space.  SSA  dictates  highly  sen¬ 
sitive  systems  track  small  objects  at  high  orbit  round  the  clock.  With  ground  based 
systems,  large  telescopes  are  often  limited  by  weather  and  atmospheric  conditions. 
The  fielding  of  large  optical  systems  with  sensitive  pointing  control  mechanisms  are 
extremely  costly.  It  would  be  beneficial  to  augment  existing  assets  with  less  expen¬ 
sive  apertures  that  provide  equivalent  or  better  resolution.  If  less  expensive  systems 
could  be  fielded,  additional  coverage  and  surveillance  capability  can  be  deployed  with 
limited  defense  budgets. 

In  addition  to  SSA,  persistent  surveillance  from  space  continues  to  be  an  emerg¬ 
ing  area  of  need.  A  broad  range  of  military  and  non-military  applications  dictate  per¬ 
sistent,  on-going  image  collection  with  large  service  times.  Atmospheric  turbulence  is 
not  a  limiting  factor  for  space-borne  sensors.  However,  aperture  size  and  service  time 
limit  capability.  Applications  ranging  from  reconnaissance  to  geological,  environmen¬ 
tal,  and  agricultural  surveys  require  persistent  coverage.  Improvement  in  satellite 
imagery  resolution  is  demanded  and  this  could  theoretically  be  achieved  at  higher  or¬ 
bits  to  include  geostationary  orbit.  It  is  proposed  that  large,  light-weight,  synthetic 
apertures  hosted  on  multiple  micro-satellites  can  bring  about  improved  resolution 
at  higher  orbits  with  less  cost  and  launch  weight.  Because  of  these  proposals,  laser 
speckle  non-optical  imaging  is  revisited  for  this  research.  The  next  section  provides 
an  technical  overview  of  speckle  formation  with  an  idealized  scenario. 

1.2  Technical  Motivation 

The  following  idealized  scenario  is  considered:  A  monochromatic  laser  uniformly 
illuminates  a  distant  object  producing  a  reflected  field  at  the  observation  plane.  The 
illuminating  field  is  assumed  to  be  either  constant  or  Gaussian  in  amplitude  and  may 
be  perturbed  and  attenuated  by  propagation  effects.  The  reflected  field  is  observed 
with  either  a  lens  system  (Fig.  1.1|  or  a  detector  array  without  a  focusing  lens 
(Fig.  11.21).  The  technique  of  detecting  the  reflected  laser  light  with  a  detector  array 
without  a  lens  has  been  referred  to  by  several  names  throughout  the  literature.  It 
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has  been  described  as  pupil  plane  imaging,  lensless  coherent  imaging,  non-imaged 
laser-speckle,  and  image  correlography.  The  last  term  refers  to  the  correlogram  (or 
object  autocorrelation)  that  can  be  formed  from  the  speckle  pattern  [44], 


Object 


Figure  1.1:  Diagram  of  Laser-Speckle  Imaging 


Laser  Source  -x( 


Pupil-Plane 
Detector  Array 


Figure  1.2:  Diagram  of  Lensless  Laser-Speckle  Imaging  [44] 


This  scenario  involves  illumination  with  a  coherent  source.  The  reflected  laser 
light,  when  observed,  produces  a  high-contrast,  granular,  “speckle”  pattern  directly 
related  to  the  roughness  of  the  object  surface.  Speckle  formation  is  observed  in  many 
applications  to  include  synthetic-aperture  radar,  ultrasound  medical  imaging,  and  co¬ 
herent,  visible  light  imaging  [19].  Figure fI751  demonstrates  an  example  speckled  image 
(A)  and  an  example  non-imaged,  laser  speckle  pattern  (B).  The  observed  speckle  in¬ 
tensity  pattern  is  different  for  observations  with  lens  and  lensless  systems.  The  laser 
speckle  observed  by  a  lens  system  is  essentially  a  noisy,  grainy  image  as  demonstrated 
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A  B 

Figure  1.3:  Example  Speckle:  [A]  Image  Formed  with  Coherent  Illumination,  [B] 
Non-imaged  Laser  Speckle  Pattern 


by  Fig.  1.3A-  However,  the  laser  speckle  pattern  produced  by  a  lensless  system  has 
no  visible  connection  to  the  object  and  just  appears  as  a  random  pattern  as  demon¬ 
strated  by  Fig.  1.3B.  This  random  pattern,  while  not  an  image  of  the  original  object, 
does  contain  information  about  the  object  embedded  in  the  statistical  nature  of  the 
speckle  pattern.  This  speckle  effect  is  not  observed  when  the  illumination  source  is 
incoherent.  With  incoherent  illumination,  speckle  is  not  observed  because  the  very 
nature  of  incoherent  light  removes  the  constructive  and  destructive  interference  ef¬ 
fects.  Often,  speckle  is  a  nuisance  and  must  be  suppressed  or  removed;  however,  this 
research  further  investigates  the  use  of  laser  speckle  patterns  for  object  recovery  in 
LADAR  remote  sensing  applications. 

In  this  research,  we  consider  the  coherent  illumination  and  free-space  propa¬ 
gation  geometry  depicted  in  Fig.  Ill  Ell  The  lensless  detector  array  is  placed  at  the 
observation  plane  (or  pupil  plane).  The  observed  speckle  pattern  occurs  when  the 
reflecting  surface  roughness  is  on  the  order  of  the  illuminating  wavelength.  We  can 
consider  the  Huy  gens- Fresnel  principle  to  explain  this  phenomenon.  Every  point 
on  the  reflecting  surface  is  a  unique  point  source  contributing  to  or  interfering  with 
the  overall  wavefront.  Due  to  the  relative  reflectivity  and  path  differences,  each  re¬ 
flecting  point  source  produces  random,  additive  phase  contributions  to  the  reflected 
wavefront.  These  random  contributions  serve  to  create  constructive  and  destructive 
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interference  leading  to  the  observed  speckle  pattern.  The  statistical  properties  of  the 
speckle  depend  upon  the  coherence  of  the  illuminating  light  and  the  random  proper¬ 
ties  of  the  surface  and  transmission  medium  [6].  Therefore,  the  observed  intensity  of 
the  reflected  held  is  directly  related  to  the  reflecting  object. 

y  v 


Figure  1.4:  Laser  Speckle  Geometry  [16] 

A  primary  technical  motivation  for  considering  the  correlography  technique  is 
the  ability  to  collect  object  information  without  the  blurring  effects  of  the  atmosphere. 
This  will  be  further  discussed  in  Section  4.2.21  Next  we  detail  the  mathematical  model 
under  consideration  for  the  idealized  scenario. 

1.2.1  Basic  Mathematical  Model.  The  electric  held  at  the  object  is  excited 
by  an  ideal,  single-mode  laser  with  a  purely  monochromatic  oscillator  with  known 
amplitude  S,  known  frequency  u,  time  t,  and  unknown  but  absolute  phase  <p.  We 
will  assume  the  laser  beam  is  linearly  polarized.  This  single-mode  laser  light  scenario, 
detailed  by  Goodman  [17],  is  modeled  as  a  random  process  that  is  both  stationary 
and  ergodic.  The  signal  in  this  ideal  scenario  is  represented  by 

u[t)  =  S  cos[2nut  —  ip],  (1.1) 

The  starting  point  with  this  simple  model  assumes  a  constant  amplitude  during  the 
duration  of  the  pulse  and  a  constant  phase  associated  with  a  sufficiently  long  coherence 
time  of  the  laser. 
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Using  scaler  theory  and  complex  notation,  we  consider  random  contributions 
from  a  large  number,  N,  of  independent  scattering  sources.  Specifically,  at  a  point 
on  the  observation  plane,  the  complex  held  is  the  result  of  a  random  phasor  sum. 
The  mathematics  of  the  random  phasor  sum  found  in  optics  is  covered  extensively 
by  Goodman  [17].  Equation  11.21  describes  the  random  phasor  sum  where  a *.  is  the 
random  amplitude  and  9 *.  is  the  random  phase  of  the  kth  phasor: 


ae ^  = 


Vn 


N 

E 

k= 1 


aye- 


jSk 


(1.2) 


The  monochromatic  illumination  produces  a  held  distribution  dependent  upon 
the  laser  beam  distribution  and  object  rehectivity.  The  rehected  object  held,  /(x), 
at  the  object  plane  is  described  with  complex  notation  as: 


/(x,t)  =  a(x,t)  •  e^(x),  (1.3) 

where  x  =  (x,y)  is  a  two  dimensional  coordinate  vector  in  the  object  plane,  a(x)  is 
the  object  held  amplitude  related  to  the  object  surface  rehectivity,  and  </>(x)  is  the 
phase  directly  related  to  the  object  surface  height  prohle.  Because  of  object  surface 
roughness,  the  rehected  phase  is  modeled  as  uniformly  distributed,  U[-vr,7r];  0(xi) 
and  0(x2)  are  independent.  This  produces  a  random  phasor  sum  at  the  detector 
plane  [17,19].  For  the  proposed  LADAR  system,  it  is  assumed  the  rehected  laser  pulse 
is  integrated  (summed)  at  the  detector  system  for  the  entire  pulse  width.  Possible 
pulse  width  changes  and  effects  due  to  oblique  reflection  angles  are  ignored. 

If  slight  changes  to  the  geometry  and  environment  for  each  laser  pulse  are  con¬ 
sidered,  the  object  held  /(x,  t)  becomes  a  random  process.  The  illuminating  laser 
may  move  slightly,  the  atmosphere  may  change  the  illuminating  laser  light  propa¬ 
gation  from  pulse  to  pulse,  and  the  object  may  move  slightly.  Slight  geometry  and 
environment  changes  produce  a  unique  rehected  held  for  each  laser  pulse.  Therefore, 
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the  reflected  field  can  be  described  in  statistical  terms.  Eqn.  1.41  demonstrates  the 
expected  value  of  the  reflected  field  is  zero,  where  E[-}  is  the  Expected  Value  operator. 


£[a(x,f)e^(x)] 
a(x,t)  -E[e^(x)] 


a(x,  t)  ■  0 


cos  0(x)  +  jsin0(x)  ]d(p 


0 


(1.4) 


The  autocorrelation  of  the  object  field  is  detailed  in  Eqn.  3  .5  where  a2(x)  =  o(x)  is  the 
incoherent  object  intensity  or  brightness  function  and  <5(x)  is  Dirac’s  delta  function. 
The  expected  value  operator  in  Eqn.llll.5lis  zero  everywhere  except  when  0(xi)  =  </>(x2) 
because  statistical  independence  in  the  phase  at  each  spatial  point  on  the  object 
surface  is  assumed.  The  reflected  field  is  described  as  a  circular  Gaussian  process 
which  provides  some  unique  properties  when  studying  the  laser  speckle  patterns  in 
the  observation  plane. 


-E[/(xi,t).f  (x2,t)] 


E[a(x1,  t)ei0(xi)a(x2,  t)e“^(x2)] 
a(xi, t)a(x2,  t)  ■  E[e’^Xl')e~:’^X2'>] 
a(xi,  t)a(x2,  t)  ■  <5(xi  —  x2) 
o(xi)  •  5(xi  -  x2) 


(1.5) 


The  reflected  field  at  the  observation  plane,  F(u,  v ),  in  the  far  field  or  Fraunhofer 
region  is  related  to  the  object  field  by  the  scaled  Fourier  Transform  [18].  This  is 
described  in  Eqn.  11.61  where  k  =  2n/\  is  the  wave  number,  A  represents  the  optical 


wavelength  of  the  monochromatic  laser  and  z  represents  the  propagation  distance 
between  the  object  and  the  observation  plane. 


F(u,v )  = 


ejkzej£(u2+v2)  roo  j-  oo 


jXz 


'  — oo  J  — oo 


2  7T 

f(x,y)exp{-j  —  (xu  +  yv)}dxdy  (1.6) 


At  optical  wavelengths  the  electric  field  cannot  be  observed  directly;  however, 
the  optical  power  or  intensity  can  be  measured  via  optical  detectors.  Following  Good¬ 
man’s  framework,  the  intensity  is  defined  as  the  squared  modulus  of  the  analytic  signal 
representation  of  the  field  [17,18].  The  intensity  at  the  observation  plane  is  described 
by 


J(u)  =  |^{/(x)}|2  (1.7) 

where  J-Xz  is  the  two-dimensional,  scaled  Fourier  Transform  detailed  in  Eqn.  11761 
J(u)  is  corrupted  by  speckle  noise  due  to  the  random  phasor  sum  produced  by  the 
object  surface  roughness.  Note  that  Eqn.  11.71  describes  the  modulus  squared  of  the 
Fourier  Transform  of  the  field  distribution  which  is  related  to  the  held  autocorrelation, 
.Rf(xo),  via  the  Fourier  Transform  due  to  the  Autocorrelation  Theorem  [18]: 

|  ^{/(x)}  |2=^{fl,(x0)}.  (1.8) 

Unfortunately,  the  complete  spectrum  cannot  be  observed  due  to  noise  and  the  finite 
extent  of  our  detection  systems.  The  observed,  speckled  intensity  data,  J0(u),  is 
modeled  as 


J0(u)  =  [/( u)  +  ra(u)]  •  A(u) 


(1.9) 
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where  n(u),  photon  or  shot  noise,  is  a  zero  mean  noise  such  that  the  observed  in¬ 
tensity  (conditioned  on  the  average  photon  values)  has  a  Poisson  distribution  with  a 
mean  equal  to  the  intensity  without  photon  noise.  This  noise  is  introduced  by  the 
random  arrival  time  of  the  photons  emerging  from  the  electric  field  onto  the  detec¬ 
tor.  The  mean  intensity,  observed  in  photon  counts,  is  itself  a  random  process  with 
an  exponential  distribution  [17].  A(u)  is  the  aperture  function  denoting  the  region 
where  the  speckle  pattern  is  physically  recorded;  71  (u)  =  1  for  the  points  within  the 
measurement  aperture  and  A(u)  =  0  elsewhere.  Without  a  traditional  lens  aperture, 
the  detector  array  is  itself  the  limiting  aperture  capturing  a  finite  portion  of  the  re¬ 
flected  field.  For  systems  designed  to  recover  an  image  of  the  original  object  from 
the  observed  laser  speckle,  it  is  important  the  laser  speckle  images  are  statistically 
independent.  Practically,  this  condition  is  easily  achieved. 

To  complete  the  description  of  the  idealized  scenario,  a  few  additional  assump¬ 
tions  and  limitations  must  be  explored.  This  non-imaged  laser  speckle  scenario  can  be 
considered  without  atmospheric  disturbance  or  where  the  atmospheric  turbulence  is 
modeled  as  a  phase  screen  directly  over  the  remote  object.  Each  pupil  plane  detector 
observes  a  beam  limited  scenario  with  the  array  producing  a  single  laser  speckle  im¬ 
age  per  pulse.  This  non-imaged  pupil  plane  intensity  is  observed  over  a  large  number 
of  laser  pulses.  The  observation  period  for  each  laser  pulse  is  assumed  to  be  within 
the  coherence  time  of  the  laser.  In  addition,  background  light  is  assumed  to  be  com¬ 
pletely  filtered  out  and  only  the  reflected  laser  light  is  observed.  The  detector  array 
is  sufficiently  in  the  far-held  for  the  Fraunhofer  approximation  to  hold.  The  pupil 
plane  array  is  considered  perpendicular  to  the  object  plane.  The  pupil  plane  and  the 
object  is  fixed  during  the  observation  period  with  sufficiently  minor  changes  pulse  to 
pulse  to  create  statistically  independent  laser  speckle  patterns. 

1.2.2  Atmospheric  Effects  on  Laser  Speckle  Imaging.  Imaging  and  optical 
propagation  through  turbulent  medium  has  been  extensively  studied  [1,35].  At¬ 
mospheric  turbulence  effects  on  laser  speckle  imaging  or  imaging  correlography  are 
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briefly  reviewed  here.  As  previously  stated,  a  single  layer  of  atmospheric  turbulence 


encountered  at  the  observation  plane  has  little  effect  on  non-imaged  laser  speckle, 
only  layered  atmosphere  need  be  considered.  Also,  a  thin  layer  of  atmosphere  at  the 
object  plane  has  little  effect  on  nonimaged,  speckle  intensity  measurements.  This  is 
explored  mathematically  along  with  an  extended,  thick  atmosphere. 


1.2. 2.1  An  Atmospheric  Layer  at  the  Observation  Plane.  If  the  sce¬ 


nario  of  a  mountaintop  observation  system  looking  skyward  at  zenith  is  considered, 
a  single,  thin  layer  of  atmosphere  at  the  pupil  plane  is  a  reasonable  and  often  used 
model  [1],  For  mathematical  and  experimental  convenience,  extensive  research  has 
been  detailed  on  how  to  analytically  replace  a  turbulent  region  with  an  equivalent 
thin,  random  screen  [35].  For  a  satellite  at  geostationary  orbit,  the  ratio  of  the  tur¬ 
bulent  layer  thickness  to  the  overall  propagation  distance  is  very  small.  For  this  case, 
the  turbulent  layer  is  considered  a  “thin”  phase  screen.  In  this  model,  only  the  phase 
of  the  optical  wave  is  distorted  by  the  turbulent  region  and  not  the  amplitude  [1]. 
Amplitude  effects  are  predominantly  observed  from  propagation  distance.  For  this 
scenario,  phase  aberration  effects  are  modeled  as  a  single  additive  phase  parameter, 
^(u).  The  phase  aberrated  field,  G(u)  at  the  observation  plane  is  described  as 


G(u)  =  F(  u)ei*<u), 


(1.10) 


where  F( u)  is  the  reflected  field  after  Fraunhofer  propagation.  The  intensity  (without 
noise)  is 


J(u)=|G(u)|2 


(1.11) 
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In  simple  terms,  Eqn.  1.111  describes  why  the  phase  aberration  due  to  atmospheric 
turbulence  at  the  observation  plane  has  no  effect  on  the  observed  intensity.  Speckle 
noise,  detector  array  effects,  system  noise,  and  scintillation  due  to  layered  atmosphere 
are  not  overcome  in  a  single  laser  speckle  observation. 

1.2. 2. 2  An  Atmospheric  Layer  at  the  Object  Plane.  Consider  the 
scenario  of  a  satellite-based  sensing  system  observing  the  Earth’s  surface.  In  this 
case,  a  single  thin  layer  of  atmosphere  at  the  object  plane  is  a  simplified,  yet  realistic 
model.  This  is  the  same  model  as  the  previous  section;  however,  the  sensing  system 
and  the  remote  target  are  reversed.  Again,  the  relatively  lengthy  propagation  distance 
above  the  atmosphere  as  compared  to  the  short  propagation  distance  through  the 
atmosphere  is  encountered.  Much  of  the  atmospheric  turbulence  effects  is  due  to  the 
dense  atmosphere  close  to  the  Earth’s  surface.  In  considering  the  Earth’s  atmosphere 
as  a  thin  phase  screen  at  the  object  plane,  the  phase  aberrated  held  is  denoted  by 
3(x)- 


g(x)  =  /(x)e^W 

=  a(x)e^(x)e^(x) 

=  a(x)eM+^x) 

=  a(x)e^'(x)  (1.12) 


From  Eqn.  11.121  the  additive  phase  abberation  only  adds  another  random  number  to 
the  random  phase  of  the  object  held,  /(x),  due  to  the  rough  object  surface.  The 
object  phase  is  a  random  variable  uniformly  distributed  between  [ — 7r,  7t]  .  Adding  a 
zero  mean,  random  number  to  a  uniformly  distributed  random  number,  {0  G  [— vr,  7r] } , 
yields  another  uniformly  distributed  number  0  G  [— 7r,7r]}  due  to  phase  wrapping. 
Therefore,  the  result  is  identical  to  the  case  without  atmosphere. 
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1.2. 2. 3  An  Extended  Atmosphere.  Rayleigh- Sommer feld  Propagation 
theory  [18]  is  used  to  address  the  more  complicated  extended  atmosphere  scenario. 
The  geometry  for  this  propagation  is  depicted  in  Fig.  11.51  The  Rayleigh- Sommer  feld 
diffraction  equation  is: 


F(u)  =  n 


/(x) 


exp 


-j'27rx(x,  u) 

A 


Ydx. 


(1.13) 


X(x,  u) 

Here,  y(x,  u)  is  the  propagation  distance  between  any  two  points  on  the  object  and 
observation  planes,  £  is  the  finite  object  plane,  spatial  area  (or  aperture  in  diffraction 
theory  terms)  and  T  is  the  obliquity  factor.  T  goes  to  unity  as  the  geometry  angles 
are  small  or  the  observing  plane  is  far  from  the  diffracting  aperture.  The  primary 
emphasis  for  exploring  Rayleigh-Sommerfeld,  Eqn.  11.131  does  not  have  the  simplifying 
assumptions  of  the  Fraunhofer  propagation. 


Object  Observation 


Figure  1.5:  Geometry  for  Propagating  Through  Turbulence 


Treating  each  point  on  the  object  plane  as  an  individual  point  source  emitting  an 
optical  wave,  the  wavefront  undergoes  a  propagation  path  length,  y  (x,  u)  dependent 
upon  spatial  position  in  both  the  object  and  observation  planes.  y(x,  u)  is  calculated 
according  to: 


x(x,y:u:v)  =  \J  z1  +  (x  —  u)2  +  (y  —  v)2. 


(1.14) 
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Note  when  z  is  large  (the  Fraunhofer  region)  and  two  simplifying  approximation  are 
made,  the  Rayleigh- Sommerfeld  diffraction  equation  simplifies  to  the  Fraunhofer  prop¬ 
agation  model  [18].  These  are  (1)  x  ~  z  for  the  term  not  found  in  the  exponent  and 
(2)  for  the  x  found  in  the  exponent,  the  binomial  expansion  to  Eqn.  1.141  is  approxi¬ 
mated  with  only  the  first  two  terms.  Returning  to  Rayleigh- Sommerfeld  propagation, 
the  atmospheric  effects  on  the  propagation  path,  y(x,  u),  must  be  included. 

The  extended  atmospheric  turbulence  produces  rapid  and  random  changes  in 
the  refractive  index  [1],  Using  the  first  order,  Rytov  approximation  (and  weak  fluc¬ 
tuation  conditions),  the  turbulence  produces  a  complex  phase  perturbation  on  the 
optical  wave  [1].  This  effectively  delays  propagation  or  produces  phase  aberrations 
for  the  propagating  wavefront.  The  phase  aberration,  ^(x,  u,  z),  essentially  changes 
the  path  length  by  the  relationship 


U\a(x,u) 
c(x,  u,  z) 


(1.15) 


where  v  is  the  frequency  of  the  light,  c  is  the  speed  of  light  for  a  given  refractive 
index,  and  xa  is  the  change  in  path  length.  Adding  the  path  length  change  due  to 
the  phase  aberration  to  Eqn.  1.131  yields  Eqn.  1.161 


F (u)  =  ^ 


/(x) 


(x(x,  u)  +  Xa(x,u)) 


exp 


-j2tt(x(x,  u)  +xa(x,u)) 

A 


dx  (1.16) 


The  path  length  change,  due  to  the  phase  aberration  caused  by  turbulence,  is  a  func¬ 
tion  of  both  object  plane  and  observation  plane  coordinates.  The  random  perturba¬ 
tions  cause  path  dependent  changes  to  the  optical  wavefront.  From  this  it  is  reasoned, 
Eqn.ll.Thlcan  not  be  simplified  to  a  scaled  Fourier  Transform.  The  layered  atmosphere 
causing  scintillation  is  a  much  more  difficult  problem.  Propagation  through  extended 
turbulence  corrupts  the  Fourier  magnitude  of  the  object  held.  Extended  atmosphere 
and  scintillation  effects  on  laser  correlography  were  explored  analytically  and  exper- 
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imentally  by  Elbaum  et  al.  [10].  The  authors  concluded  under  certain  conditions 
the  random  apodization  effect  will  be  negligibly  small  for  a  ground-based  receiver 
observing  a  near-zenith  exoatmospheric  target.  This  research  will  only  consider  the 
simple,  atmospheric  model  of  a  thin  phase  screen  either  at  the  object  or  observation 
planes.  Although  previous  researchers  in  this  area  acknowledge  effects  due  to  high 
altitude  turbulence  and  scintillation,  a  thin  atmosphere  or  screen  is  often  assumed 
for  simplification  or  implied  with  the  data  model  [11,23,25,39,44],  The  thin  screen 
assumption  is  continued  throughout  this  dissertation. 

1.3  Research  Contributions 

This  dissertation  provides  three  primary  research  contributions: 

•  Statistical  Analysis  of  Correlography  Data.  The  observed  laser  speckle  is  trans¬ 
formed  via  a  post-processing  technique  that  produces  noisy  autocorrelations. 
The  randomness  of  the  transformed  data  can  be  closely  modeled  by  the  nega¬ 
tive  exponential  distribution. 

•  Two  Iterative  Algorithms  Using  Polarimetric  Data.  Iterative  solutions  for  max¬ 
imum  likelihood  and  maximum  a  posteriori  estimators  are  developed  using  the 
Expectation  Maximization  technique. 

•  A  Theoretical  Bound  on  Resolution  for  Correlography.  A  Cramer-Rao  Lower 
Bound  for  resolving  two  closely  spaced  point  sources  is  presented  using  the 
negative  exponential  noise  model  and  polarimetric  data  models. 

1-4  Dissertation  Overview 

This  document  is  divided  into  seven  chapters  and  contains  three  appendices. 
Chapter  2  presents  relevant  technical  background  information  as  the  theoretical  basis 
for  this  research.  Previous  research  in  phase  retrieval  and  non-imaging  laser  speckle 
is  highlighted. 
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Chapter  3  provides  the  relevant  mathematical  models  and  theory  that  form  a 
foundation  for  this  research.  Although  most  of  this  material  may  be  found  elsewhere, 
it  provides  foundational  elements  behind  the  later  chapters  to  include  proposed  sensor 
models  and  statistical  framework. 

Chapter  4  provides  new  theoretical  work  associated  with  computational  algo¬ 
rithms  for  solving  the  phase  retrieval  problem  using  polarization  diversity.  Two  po- 
larimetric  algorithms  are  presented  differentiated  by  data  collection  schemes. 

Chapter  5  details  a  theoretical  resolution  bound  for  correlography  using  polar¬ 
ization  diversity.  The  bound  itself  is  not  predictive  of  specific  algorithm  performance; 
however,  it  demonstrates  theoretical  improvement  provided  by  adding  polarization 
diversity.  A  Cramer-Rao  Lower  Bound  is  presented  for  both  unpolarized  and  po¬ 
larized  data  scenarios.  The  comparison  of  theoretical  bounds  demonstrates  relative 
improvement  provided  by  polarization  diversity. 

Chapter  6  presents  results  and  analysis  from  computer  simulation  and  labora¬ 
tory  experiments.  A  representative  subset  of  simulation  and  experimental  results  are 
presented  to  support  key  research  findings  and  contributions.  This  chapter  provides 
comparative  results  between  polarimetric  and  non-polarimetric  algorithms. 

Chapter  7  concludes  the  main  document  by  providing  an  overall  summary  of 
research  activities,  a  summary  of  key  findings,  and  recommendations  for  subsequent 
research.  This  is  followed  by  the  appendices  that  provide  some  of  the  detailed  math¬ 
ematical  proofs. 
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II.  Background 

This  chapter  describes  the  theoretical  background  and  previous  research  work  for  ob¬ 
ject  recovery  from  laser  speckle  images.  The  ultimate  goal  of  object  recovery  is  the 
synthesis  of  the  best  possible  object  image  from  a  limited  data  set.  With  intensity 
measurements,  the  phase  of  the  original  object  field  is  unknown.  In  fact,  the  phase 
information  is  lost  during  the  detection  process,  presenting  an  ill-posed  problem.  In 
the  LADAR  application  considered  in  this  research,  an  estimate  of  the  object  auto¬ 
correlation  is  recorded.  The  Autocorrelation  Theorem  relates  the  autocorrelation  to 
the  object’s  Fourier  Magnitude  via  the  Fourier  Transform  [18].  Unfortunately,  com¬ 
plete  information  about  the  original  object  is  not  observed.  However,  with  iterative 
computing  techniques,  it  has  been  shown  that  an  image  solution  relatively  similar  the 
original  may  be  produced.  This  technique  is  often  described  as  phase  retrieval  and  is 
employed  in  several  applications  found  in  image  and  signal  processing.  The  goal  of 
this  research  is  to  provide  improvement  for  this  specific  phase  retrieval  problem  with 
laser  speckle  images. 

2.1  Averaging  Laser  Speckle  Patterns 

Often  with  signals  corrupted  by  random  noise,  the  first  course  of  action  is  to 
average  a  large  number  of  independent  signal  realizations  in  hopes  of  minimizing  the 
noise  effects.  In  fact,  with  laser  speckled  images  (formed  with  a  lens),  averaging  a 
large  number  of  registered,  independent,  speckled  images  does  reduce  speckle  effects. 
There  are  several  techniques  for  speckle  suppression  in  optical  imaging  [19].  However, 
averaging  nonimaged  laser  speckle  patterns  does  not  yield  the  same  result.  The  en¬ 
semble  average  of  independent,  nonimaged  laser  speckle  intensity  patterns  converges 
to  a  constant,  C, 


K 

lim  K  1  \  IQk( u)  =  U,  (a  constant), 

iC— >o o  ^ ' 

k=  1 


(2.1) 
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where  K  is  the  number  of  observed  frames  and  IQk  is  the  kth  frame  of  observed, 
speckled  intensity  data.  This  is  proven  in  Appendix  |AJ  Because  of  this  result,  we 
must  look  to  another  post  processing  method  to  recover  the  object. 

2.2  Object  Autocorrelation  from  Intensity  Data 

Due  to  speckle,  photon  noise,  and  detector  array  effects,  a  single  laser  speckle 
image  is  not  useful  for  object  recovery.  The  random  speckle  pattern  contains  statis¬ 
tical  information  about  the  object  surface,  though  insufficient  for  human  perception. 
However,  Idell  et  al.  pointed  out  the  ensemble  average  of  the  magnitude  squared  of 
the  Fourier  transform  of  the  intensity  data  converges  to  a  result  directly  related  to 
the  object  autocorrelation  [25].  We  will  call  this  operation  the  Idell  function  and 
describe  the  results  in  Eqn.  |2.2[ 

K 

lim  l^"1{4fe(u)}|2  ]  =b\h(x)\2 +  c[R0(x)]*\h(x.)\2  (2.2) 

K —XX)  1 

k= 1 

where  b  and  c  are  constants,  K  is  the  number  of  independent  speckle  realizations, 
i?Q(x)  is  the  autocorrelation  of  the  object  intensity,  |h(x)|2  is  the  incoherent,  intensity 
impulse  response  of  the  detector  array  and  *  represents  convolution  operation.  The 
object  autocorrelation  is  blurred  by  the  impulse  response  of  the  detector  array  in  a 
similar  fashion  as  an  optical  system  blurs  an  image  via  diffraction  effects.  Because  of 
this,  the  resulting  autocorrelation  is  considered  diffraction-limited  [13].  The  impulse 
response  of  the  detector  array  should  be  known  providing  the  ability  to  produce  a 
very  good  estimate  of  the  object  autocorrelation.  The  result  in  Eqn.  12.21  ignores 
photon  noise  added  during  the  detection  process.  The  proof  of  this  result  is  found  in 
Appendix  [B] 

Through  the  Wiener-Khinchin  Theorem,  the  object  autocorrelation  is  related 
to  the  power  spectral  density  via  a  Fourier  Transform  relationship.  Here,  the  phase 
retrieval  problem  is  presented  as  only  the  Fourier  magnitude  is  estimated.  In  the 
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next  section  of  this  chapter,  methods  for  recovering  the  Fourier  phase  of  the  unknown 
object  are  discussed. 

2.3  Phase  Retrieval 

The  problem  of  phase  retrieval  comes  about  when  measured  data  presents  only 
the  magnitude  of  the  signal’s  Fourier  Transform.  In  the  LADAR  application  presented 
in  this  paper,  the  intensity  or  magnitude  squared  is  collected.  Without  constraints  or 
knowledge  of  the  signal,  the  loss  of  the  phase  is  irreversible  [21].  However,  from  the 
body  of  research  available  on  phase  retrieval  problems,  with  certain  constraints  and 
a  priori  knowledge  of  the  signal,  a  recovery  is  possible. 

The  idea  of  synthesizing  images  from  non-imaged  laser-speckle  is  not  new.  In 
1987,  Idell  et  al.  proposed  a  method  for  recovering  unspeckled  images  and  demon¬ 
strated  this  technique  possible  with  computer  simulation  [25].  It  is  this  seminal  paper 
where  the  starting  point  of  our  research  is  found.  An  estimate  of  the  object  auto¬ 
correlation  is  formed  from  many  laser  speckle  images.  From  this,  a  phase  retrieval 
algorithm  is  employed  to  produce  an  estimate  of  the  object.  Several  phase  retrieval  al¬ 
gorithms  have  been  proposed  in  published  literature  that  support  our  LADAR  remote 
sensing  application.  A  review  of  these  algorithms  is  presented  next. 

2-4  Phase  Retrieval  Algorithms 

Dainty  and  Fienup  detail  a  thorough  review  of  several  phase  retrieval  methods 
to  include  Newton- Raphson,  gradient  search,  and  iterative  methods  [7,12],  The  phase 
retrieval  problem  presents  itself  in  many  applications  to  include  optical  images  formed 
with  incoherent  illumination.  For  the  application  considered  here,  a  short  review  of 
applicable  phase  retrieval  algorithms  is  presented. 

2-4-1  Gerchberg- Saxton  Algorithm.  Gerchberg  and  Saxton  suggested  iter¬ 
ating  between  two  sets  of  data  that  are  related  by  the  Fourier  Transform  [15].  They 
suggest  simultaneously  recording  intensity  measurements  in  both  the  image  and  pupil 
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planes.  The  algorithm  starts  with  an  initial  phase  guess  generated  by  a  uniform  ran¬ 
dom  number  generator  (—7 r  to  tt).  The  guessed  phase  is  multiplied  to  the  measured 
amplitude  in  the  image  plane  and  the  Fourier  Transform  is  taken  mapping  the  data 
to  the  pupil  plane.  The  computed  phase  from  this  operation  is  multiplied  to  the 
amplitude  of  the  pupil  plane  data.  Again,  a  Fourier  transform  is  performed  mapping 
the  manipulated  data  set  back  to  the  image  plane.  The  computed  phase  in  this  op¬ 
eration  is  applied  to  the  amplitude  of  the  image  plane  data.  This  iterative  operation 
continues  until  operator  intervention  or  the  desired  level  of  error  is  reached.  Gerch- 
berg  and  Saxton  demonstrated  the  estimate  error  decreases  or  remains  constant  for 
each  iteration.  This  is  an  attractive  property  of  the  algorithm;  it  avoids  a  diverging 
solution. 

2-4-2  Fienup’s  Algorithms.  Fienup  details  two  phase  retrieval  algorithms 
descending  from  the  work  of  Gerchberg  and  Saxton.  The  first,  Fienup’s  Error  Reduc¬ 
tion  method  is  essentially  a  generalized  form  of  the  Gerchberg-Saxton  phase  retrieval 
algorithm  [11, 12]  and  is  often  quoted  in  literature.  The  algorithm  essentially  iter¬ 
ates  between  object  and  Fourier  domains  where  known  constraints  are  applied  to  the 
data  before  continuing  with  the  next  iteration.  For  the  object  domain,  the  object 
is  assumed  to  be  positive  and  within  a  known  observation  region,  called  object  sup¬ 
port.  In  the  beam-limited,  LADAR  application,  the  support  region  is  provided  by 
the  extent  of  the  illuminating  laser  beam.  For  the  Fourier  domain,  the  modulus  of 
the  Fourier  Transform  of  the  object  field  is  known  from  the  observed  data.  For  the 
LADAR  application  considered  here,  this  results  from  the  square  root  of  the  Fourier 
Transform  of  the  noisy  (estimated)  autocorrelation  data.  Figure  2.1  depicts  a  block 
diagram  of  Fienup’s  Error  Reduction  algorithm  [11], 

An  initial  guess  of  the  object  is  required  for  the  algorithm.  Without  any  knowl¬ 
edge  of  the  object,  the  guess  can  be  from  an  infinite  number  of  choices.  However,  the 
extent  of  the  object  or  support  region  may  be  known  or  assumed.  Often,  the  initial 
guess  is  chosen  to  be  the  support  region  itself  perturbed  with  random  noise.  In  a 
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Figure  2.1:  Block  Diagram  of  Fienup’s  Error  Reduction  Method  [11] 

later  paper  [12],  Fienup  proposes  stopping  criteria  to  be  the  squared  error  between 
the  Fourier  constraint  and  the  computed  Fourier  Transform  (similar  in  the  object  do¬ 
main  for  successive  iterations).  This  algorithm  is  a  single  frame  algorithm;  it  operates 
on  a  single  autocorrelation  estimate.  Without  stopping  criteria,  the  algorithm  runs 
until  operator  intervention  or  the  number  of  chosen  iterations  has  been  exceeded. 

Fienup  also  presents  the  input-output  algorithm  to  speed  convergence  [11,12], 
Fienup  adapts  the  error  reduction  algorithm  to  produce  a  non-linear  approach  where 
the  input  does  not  necessarily  satisfy  the  object  domain  constraints.  Also,  the  input 
to  the  algorithm  at  each  iteration  is  not  considered  the  current  best  estimate  and 
can  be  modified.  Fienup  suggests  multiple  methods  for  selecting  a  new  estimate  of 
the  object  input  and  recommends  periodically  changing  the  selection  method  after 
several  iterations.  Fienup  also  suggests  a  hybrid  approach  where  the  error  reduction 
algorithm  is  combined  with  the  input-output  algorithm.  Although  Fienup  reports  this 
as  a  very  powerful  approach,  it  appears  ad  hoc  through  arbitrary  user  intervention. 

2-4-3  Schulz-Snyder  Algorithm.  Schulz  and  Snyder  present  a  unique  image 
recovery  algorithm  that  operates  on  nth  order  correlations  [38].  A  phase  retrieval 
algorithm  is  presented  with  n  =  2.  Schulz  and  Snyder  choose  to  maximize  a  log- 
likelihood  function  of  the  data  where  a  Poisson  model  is  selected.  Schulz  and  Snyder 
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acknowledge  that  few  applications  present  autocorrelation  data  corrupted  by  Poisson 
noise;  however,  he  states  this  model  does  enforce  positivity  and  is  similar  to  cases 
where  the  noise  is  signal  dependent.  For  applications  with  signal  dependent  noise  but 
unknown  distribution  this  algorithm  may  be  a  good  choice.  For  the  application  con¬ 
sidered  in  this  paper,  the  autocorrelation  does  not  exhibit  Poisson  noise  as  discussed 
in  Chapter  Hill  Equation  {2.31  depicts  the  Schulz-Snyder  iterative  algorithm: 


Ok+ i(x) 


Ok{x) 

2ExOk(x)  ^ 


ok(x  +  y)  +  ok(x  -  y ) 


Ro(y) 
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(2.3) 


where  ok+ 1  is  the  new  estimate  of  the  object,  ok  is  the  old  estimate  of  the  object 
from  the  previous  iteration,  R0{y)  is  the  measured  autocorrelation  (observed  data), 
and  R0{y)  is  the  autocorrelation  formed  from  the  old  object  estimate.  Using  notation 
found  within  this  document,  the  Schulz-Snyder  algorithm  is  re-written  to  be 
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where  *  is  correlation,  *  is  convolution,  onew  is  the  new  object  estimate,  oold  is  the 
object  estimate  from  a  previous  iteration,  Ra  is  the  measured  autocorrelation,  i?°ld 
is  the  autocorrelation  formed  from  the  old  estimate  of  the  object,  and  S°ld  is  the 
two-dimensional  sum  of  the  old  object  estimate. 

The  Schulz-Snyder  algorithm  also  requires  a  support  region  where  the  object  is 
known  to  exist  and  an  initial  guess.  Schulz  and  Snyder  caution  against  using  a  smooth, 
uniform  image  as  the  initial  guess  and  suggest  an  asymmetric  starting  image  [38]. 
The  initial  guess  for  the  Schulz-Snyder  algorithm  can  also  be  a  gross  starting  point 
such  as  the  support  region  perturbed  with  independent  random  variables  distributed 
uniformly  over  the  interval  [0.99,1.01].  The  Schulz-Snyder  algorithm  is  also  a  single 
frame  algorithm  and  no  stopping  criteria  is  provided.  The  Schulz-Snyder  iterative 
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algorithm  can  be  characterized  as  a  gradient  search  method  with  object  constraints 
used  to  aid  in  convergence.  The  initial  guess  (positive  with  support)  enforces  a 
solution  only  within  the  support  region  as  well  as  object  positivity. 

2-4-4  Phillips’  Algorithm.  Phillips  presents  a  statistical  approach  to  the 
phase  retrieval  problem  claiming  Gaussian  statistics  when  the  number  of  independent 
laser  speckle  patterns,  K ,  is  large  [32],  While  the  Central  Limit  Theorem  suggests  this 
is  true  under  the  right  conditions,  Chapter  [TTl]  suggests  this  is  a  poor  approximation 
for  our  LADAR  application.  From  the  Gaussian  assumption,  the  log-likelihood  is: 

L(°)  =  ~2^y)  S  -  °(x)°(x  +  y)  j  >  (2-5) 

where  o{x)  is  the  true  object,  d(y)  is  the  average  speckled  autocorrelation,  and  cr2(y) 
is  the  variance.  An  iterative  maximization  approach  is  accomplished  with  the  kth 
frame  data  model  detailed  in  Eqn.  12.61  Averaging  the  data  over  K  number  of  frames 
is  assumed  to  be  Gaussian  with  a  mean  equal  to  the  Idell  function  (Eqn.  12.2]).  A 
sample  variance  is  computed  from  the  K  observed  frames. 


dk(  y)  = 


T  1  +  n‘(u)} 


(2.6) 


Phillips  suggests  a  low-resolution  image  as  the  starting  point.  The  author  found 
the  Phillips  algorithm  to  converge  with  a  random  guess  as  the  starting  image  identical 
to  Fienup  and  Schulz- Snyder.  The  Phillips  algorithm  may  be  implemented  as  a 
multi-frame  algorithm  iteratively  operating  on  each  noisy  autocorrelation  realization 
obtained  after  each  laser  pulse.  Lastly,  Phillips  suggests  a  statistical  based  stopping 
criteria. 
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2.5  Additional  Image  Recovery  Efforts 

In  addition  to  the  above  phase  retrieval  algorithms,  additional  methods  for 
object  recovery  have  been  suggested  for  identical  or  similar  applications. 

2.5.1  Phillips  Joint  Algorithm.  Originally  proposed  by  Cain  [3],  Phillips 
and  Cain  describe  a  joint  algorithm  using  both  image  and  pupil  plane  data  together 
to  aid  image  recovery  [33].  The  application  considered  is  identical  to  the  application 
in  this  research:  a  coherent  LADAR  scenario  with  both  image  and  pupil  plane  data 
available.  The  Phillips  and  Cain  technique  essentially  combines  blind  deconvolution 
of  laser  imaging  and  phase  retrieval  of  laser  speckle  pattern  data  to  produce  the 
reconstructed  image.  A  Bayesian  method  is  presented  that  assumes  a  statistical  model 
for  the  image  and  pupil  plane  data  sets.  A  combined  joint  probability  density  function 
is  developed  and  subsequently  a  joint  log  likelihood  function.  The  log  likelihood 
function  is  then  maximized  through  an  iterative  technique  as  a  solution  for  the  object 
estimate.  Two  primary  assumptions  are  taken  to  develop  this  algorithm.  First, 
the  image  and  pupil  plane  data  sets,  though  collected  simultaneously,  are  statistically 
independent.  Second,  the  average  speckle  autocorrelation  (from  the  pupil  plane  data) 
is  approximated  as  Gaussian.  While  the  first  assumption  is  presented  without  proof, 
it  is  reasonable  if  the  image  and  pupil  plane  data  are  collected  at  angular  or  range 
offsets.  With  a  slight  change  in  angle,  the  random  phase  contributions  from  the  rough 
surface  produce  a  different  speckle  realization  at  the  image  and  pupil  plane  collections. 
The  second  assumption  is  provided  from  the  Central  Limit  Theorem  as  a  result  of 
averaging  a  large  number  of  laser  speckle  autocorrelations  in  the  processed  data.  This 
assumption  is  questioned  for  this  research  and  discussed  further  in  Chapter  ITTT1 

The  Phillips  Joint  algorithm  provides  an  important  method  as  many  remote 
sensing  applications  use  an  optical  imaging  system  or  have  an  imaging  capability 
available.  However,  the  image  recovery  result  depends  upon  the  resolution  of  the 
imaging  system  and  the  number  of  laser  speckle  patterns  collected.  The  better  reso¬ 
lution  of  the  imaging  system  and  the  more  laser  speckle  patterns  collected,  the  better 
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this  algorithm  performs.  Often,  optical  resolution  and  data  collection  time  are  pre¬ 
mium  quantities.  This  research  hopes  to  decouple  the  need  for  the  optical  imaging 
system  and  reduce  the  number  of  laser  speckle  patterns  required  through  better  sta¬ 
tistical  models. 

2.5.2  Sparse  Arrays  of  Detectors.  Fienup  and  Idell  proposed  nonimaged 
laser  speckle  pattern  collection  with  large,  yet  sparse  or  partially  filled  detection  ar¬ 
rays  [14].  This  research  extends  the  previous  work  by  the  same  authors  to  a  large 
array  developed  from  a  synthetic  array  of  subapertures.  This  partially  filled  array 
does  distort  the  recovered  image  via  a  transfer  function.  However,  Fienup  and  Idell 
demonstrate  with  sufficient  post-processing  the  transfer  function  effects  can  be  largely 
removed.  This  result  is  important  for  applications  where  a  very  large  aperture  is  re¬ 
quired;  however,  is  synthesized  with  smaller  subaperture  hosted  on  multiple  vehicles 
or  geographically  dispersed.  It  is  important  to  note  the  synthetic  array  is  not  devel¬ 
oped  by  translating  a  single  subaperture  in  time  as  in  Synthetic  Aperture  RADAR. 
For  the  application  considered  in  this  research,  multiple  subapertures  collect  data  at 
the  same  time. 

2.5.3  Other  Novel  Ideas.  Additional  ideas  are  presented  throughout  litera¬ 
ture  and  are  too  numerous  to  summarize  in  this  document.  Many  of  the  published 
phase  retrieval  algorithms  and  methods  do  not  apply  to  the  application  considered  in 
this  research.  However,  a  few  applicable  techniques  are  worthy  of  mention. 

Guizar-Sicairos  and  Fienup  address  the  effects  of  a  finite  sized  detector  array 
on  the  iterative  reconstruction  algorithms  using  Fast  Fourier  Transforms  (FFTs)  [20]. 
In  any  physical  system  implemented  to  measure  backscatter  energy  from  coherent 
illumination,  the  reflected  field  extends  beyond  the  width  of  the  detector  array.  The 
authors  suggest  this  discontinuity  in  the  truncated  intensity  data  causes  aliasing  or 
non-physical  effects  during  computer  simulation.  Also,  a  finite  object  will  have  an 
infinite  spectrum  causing  computational  problems  with  iterative  techniques.  The 
authors  propose  a  weighted  projection  in  the  Fourier  domain  to  compensate.  Guizar- 
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Sicairos  and  Fienup  state  this  technique  improves  image  quality  without  sacrificing 
resolution  and  is  robust  in  the  presence  of  noise.  For  the  research  considered  in  this 
paper,  the  proposed  physical  system  will  involve  a  beam-limited  scenario.  The  laser 
beam  illuminating  any  real  object  will  normally  have  significant  roll-off  in  power  at 
the  edges  of  the  beam.  Therefore,  a  “hard  object  edge”  from  the  object  or  assumed 
support  region  as  described  by  Guizar-Sicairos  and  Fienup  is  not  encountered.  The 
finite  support  in  our  research  is  considered  to  be  the  extent  of  the  illuminating  beam 
without  further  assumption  of  the  object  extent.  This  research  will  ensure  aliasing 
due  to  computation  is  eliminated  through  proper  sampling  in  the  observation  plane. 

In  another  research  effort,  Seldin  and  Fienup  showed  the  use  of  the  Ayers-Dainty 
two-dimensional  Blind  Deconvolution  Algorithm  applied  to  phase  retrieval  [40].  The 
Ayers-Dainty  algorithm  is  an  important  result  with  broad  implications  throughout  re¬ 
mote  sensing  and  image  processing.  Blind  deconvolution  is  attempted  when  one  seeks 
to  solve  for  two  unknown  functions  from  a  single  noisy  measurement  (see  Eqn.  |2.7j). 
This  is  often  the  case  in  imaging  where  the  point  spread  function,  g(x),  of  the  system 
is  unknown  due  to  atmospheric  distortion.  In  the  blind  deconvolution  problem  for 
imaging,  one  hopes  to  recover  the  image,  /(x),  through  iterative  computation  where 
the  point  spread  function  must  be  estimated. 

d(x)  =  /(x)  *  #(x)  +  ™(x)  (2.7) 


Seldin  and  Fienup  point  out  that  phase  retrieval  is  a  special  case  of  blind  deconvolu¬ 
tion  when  n(x)  =  0  and  g(x)  =  /*(— x)  because  of  the  relationship  of  the  autocorre¬ 
lation,  i?(x),  to  the  modulus  squared  of  the  Fourier  Transform. 
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^(x)  =  /(x)  *  /*  (~x) 

=  rH  F(u)F*( u)} 

=  T~l{  |F(u)|2  }  (2.8) 


Note  that  autocorrelation  is  a  special  case  of  the  convolution  theorem  when  /(x)  is 
convolved  with  /*(— x)  [18].  Seldin  and  Fienup  applied  the  Ayers-Dainty  algorithm 
and  assumed,  /  and  n  are  independent,  zero-mean,  Gaussian  random  processes  for  the 
Wiener  filter.  The  authors  acknowledge  most  images  do  not  satisfy  this  assumption. 
However,  they  show  the  filter  is  effective  with  the  end  result  equivalent  to  the  Fienup 
error  reduction  algorithm  in  zero  or  low  noise  cases.  While  this  approach  is  not 
a  significant  improvement  over  other  techniques,  it  does  further  demonstrate  the 
mathematical  importance  of  the  phase  retrieval  problem.  Correct  statistical  models 
are  significant  for  improving  upon  existing  phase  retrieval  algorithms. 

This  research  investigated  polarization  diversity  for  improving  the  phase  re¬ 
trieval  problem.  Multi-channel  diversity  is  not  a  new  concept.  Holmes  et  al.  applied 
several  iterative  methods  to  intensity  data  formed  from  two  illuminating  sources  sep¬ 
arated  in  frequency  [23,24].  For  this  specific  detection  scheme,  the  Expectation- 
Maximization  algorithm  approach  is  further  motivated  and  detailed  mathematically 
by  Schulz  and  Voelz  [39].  Essentially,  two  fields,  E\  and  E2,  separated  by  frequency 
illuminate  the  distant  object  and  the  fields  are  reflected  to  an  array  of  detectors  such 
that  the  individual  field  component  magnitudes  and  a  field  cross  product  is  measured. 
The  measured  intensity  in  the  pupil  plane  at  time  t  is 

J(x,t)  =  |i?i(x)  exp(jwit)  +  E2(x)  exp(jw2t)\2,  (2.9) 
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where  W\  and  w2  are  slightly  different  frequencies.  Holmes  et  al.  state  the  speckle 
held  intensities  |£i(x)]2  and  \E2(x)\2  must  be  measured  separately  in  order  to  solve 
four  equations  for  the  two  unknown  fields  E\  and  E2  from  noisy  data  A,  I2l  Ir,  Ii 
related  by: 


Ji(x)  =|Ei(x)|2  +  rn(x), 

(2.10) 

J2(x)  =\E2(x)\2  +  n2(x), 

(2.11) 

Jr(x)  =Re[  £’i(x)E,2(x)  ]  +ffi-(x), 

(2.12) 

A(x)  =Im[  E!(x)E;(x)  ]  +7ii(x), 

(2.13) 

where  ni,  n2,  nr,  7ij  represent  noise  and  *  represents  conjugation.  Holmes  et  al.  view 
this  as  a  phase  retrieval  problem  solving  for  both  reflected  fields  with  a  cross-phase 
constraint  provided  by  -Ei(x)i?2  (x).  I11  companion  papers,  the  authors  review  several 
algorithms  for  this  problem  in  both  the  root-matching  method  and  iterative  algo¬ 
rithm  method  [23,24].  The  two  separate  held  intensities,  |i?i(x)|J  and  |i?2(x)|~  are 
unrelated;  they  are  two  distinct  measurements  of  the  target.  Only  with  the  difficultly 
obtained  held  cross  product  can  this  problem  be  attempted.  The  simulated  results 
are  attractive;  however,  the  practical  implementation  with  optical  hardware  appears 
to  be  extremely  difficult  as  the  held  cross  product  must  be  isolated  from  the  com¬ 
bined  intensity  measurement.  The  held  cross  product  is  easily  obtained  at  RADAR 
frequencies  but  is  problematic  at  optical  frequencies  requiring  complex  interferomet¬ 
ric  hardware.  Obtaining  the  held  cross  product,  Ei(-x)E%(x)  introduces  an  additional 
processing  and  estimation  problem  avoided  in  this  research. 

For  the  Holmes  et  al.  two  held  mixing  approach  discussed  above,  Schulz  and 
Voelz  detail  the  theory  and  algorithm  for  the  generalized  expectation-maximization 
(EM)  method.  A  Poisson  probability  mass  function  for  the  data  is  assumed  similar 
to  the  Schulz-Snyder  phase  retrieval  algorithm  from  cross-correlations.  The  Schulz 


and  Voelz  result  is  of  interest  to  this  research  as  an  EM  algorithm  is  presented  in  this 
dissertation.  This  is  further  discussed  in  Chapters  HI  and  JVI 
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III.  Research  Foundations 


This  chapter  provides  the  sensor  and  mathematical  models  used  for  this  research  as 
well  as  some  preliminary  developments  for  the  task  of  object  recovery.  In  addition  to 
the  basic  mathematical  model  of  the  observed  and  processed  data,  a  statistical  frame¬ 
work  is  presented  as  a  basis  for  the  new  research.  A  case  is  made  for  the  exponential 
probability  density  function  describing  the  random  nature  of  the  processed  data.  This 
probability  description  is  used  throughout  the  remaining  document,  supporting  ob¬ 
ject  recovery  and  a  theoretical  resolution  bound.  Lastly,  maximum  likelihood  (ML) 
and  EM  algorithms  are  explored  using  the  exponential  distribution. 

3.1  Sensor  Model 

All  of  the  optical  detection  technology  required  to  perform  the  LADAR  task 
described  in  Chapters  U  and  II  exist  as  commercially  available  hardware  components. 
With  proper  engineering  and  system  integration,  a  cost  effective  correlography  and 
phase  retrieval  sensor  system  can  be  built  with  commercially  available  hardware  com¬ 
ponents.  Laser  engineering  and  power  considerations  must  be  included  based  on  the 
target  size  and  distance.  This  section  provides  the  proposed  sensor  model  used  to 
develop  this  research.  Individual  hardware  component  contributions  to  noise  as  well 
as  calibration  issues  are  not  considered.  Adding  polarization  sensitivity  to  the  corrcl- 
ography  system  requires  additional  considerations  detailed  below. 

3.1.1  Proposed  Sensor  Hardware.  A  large  detection  array  may  be  “synthe¬ 
sized”  using  many  individual  detection  elements  without  a  large  focusing  lens.  The 
array  may  be  built  to  be  quite  large  without  the  limitation  of  a  single  mounting 
frame  but  only  limited  by  pointing  and  integration  requirements.  The  individual  de¬ 
tector  elements  must  be  pointed  in  the  same  direction  to  eliminate  non-linear  array 
effects.  A  square,  uniform  detection  array  is  considered  for  simplicity.  With  the  use 
of  commercial  off  the  shelf  hardware  and  the  avoidance  of  a  large,  monolithic  optical 
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lens  and  refractive  mirror,  system  developers  can  obtain  significant  cost  savings  over 
traditional  optical  systems. 

A  proposed  scheme  consists  of  employing  polarizing  beam  splitters  (PBSs)  to 
produce  two  orthogonally  polarized  channels  of  observed  data.  A  single  detection 
element  consists  of  a  PBS  followed  by  a  small  lenslet  to  focus  the  optical  waveform 
and  an  optical  detection  element  such  as  an  avalanche  photodiode  (APD).  Although 
not  a  required  element,  the  lenslet  is  extremely  small  and  cost  efficient  for  focusing 
the  optical  energy  onto  the  detection  element.  The  lenslet  does  not  aid  in  forming 
a  traditional  image.  The  APD  would  be  followed  by  an  analog-to-digital  (A-D)  con¬ 
verter.  The  digital  data  would  be  stored  and  processed  using  a  computer  processor. 
Figure  15711  depicts  a  single  hardware  sub-element  of  the  proposed  array  hardware.  In 
this  manner,  a  single  array  can  be  built  for  detecting  two  orthogonal  channels. 


Figure  3.1:  Single  Hardware  Sub-element  for  Polarimetric  Correlography  System 

The  PBS  produces  two  orthogonally  polarized  data  channels,  normally  referred 
to  as  S  and  P  polarized  channels.  The  S  and  P  channels  are  attenuated  by  the 
effect  of  the  polarization  and  degree  of  transmission  through  the  PBS.  Commercially 
available  PBSs  may  be  obtained  with  transmission  efficiencies  greater  than  90  percent. 
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Computer  simulations  presented  in  Chapter  EE  assume  transmission  rates  of  100 
percent  without  significant  departure  from  true  hardware  performance. 

The  two  orthogonally  polarized  channels  allow  for  unpolarized  data  collection 
as  well  (e.g.  S  +  P);  however,  statistical  independence  between  the  polarized  and 
unpolarized  data  sets  is  a  consideration  for  system  designers.  Normally,  S  +  P  will 
not  be  statistically  independent  from  either  S  or  P  channels;  therefore,  an  alternative 
design  must  be  considered. 


Figure  3.2:  Alternative  Sub-element  for  Polarimetric  Correlography  System 

Use  of  a  traditional  beam  splitter  followed  by  polarizing  film  (or  polarization 
analyzer)  for  the  polarized  channel  may  be  employed  but  at  a  loss  of  light  levels  due 
to  the  nature  of  the  beam  splitter.  A  traditional,  non-polarizing  beam  splitter  will 
create  two  channels  but  only  half  the  original  light  is  transmitted  to  each  channel. 
The  polarized  channel  is  further  attenuated  by  the  effect  of  the  analyzer.  Figure  13.21 
demonstrates  this  system  design  using  a  traditional  beam  splitter.  It  is  assumed  here 
the  PBS  or  non-polarizing  beam  splitter  provides  only  linear  polarization  effects  and 
transmission  efficiencies  are  identical  for  the  two  output  channels. 
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3.1.2  Object  Illumination.  In  order  to  obtain  an  estimate  of  the  object 
autocorrelation  (or  equivalently  an  estimate  of  the  object’s  Fourier  magnitude),  the 
correlography  technique  requires  the  object  to  be  coherently  illuminated.  In  the 
application  considered  here,  the  object  is  either  spatially  limited  (e.g.  satellite  in 
orbit)  or  the  the  illuminating  beam  itself  provides  a  spatial  limit.  Therefore,  the 
illuminating  beam  must  be  spatially  coherent  across  the  entire  extent  of  the  beam  or 
object  surface.  This  is  achievable  with  currently  available  technology.  A  reasonable 
scenario  is  considered  with  the  geometry  shown  in  Fig.  3.3i  A  satellite  with  d,2  =  50 
meters  extent  is  orbiting  at  z  =  36,  000  kilometers  and  the  illuminating  source  has  a 
finite  extent  of  less  than  d\  =  1  meter.  The  light  originating  at  the  edges  of  the  source 
undergo  different  time  delays,  T\jc  and  r^/c.  If  the  time  difference  (r2  —  ri)/c  is  much 
less  than  the  coherence  time,  rc,  of  the  source,  then  the  spatial  coherence  (in  the  object 
plane)  is  observed  [17].  With  the  proposed  scenario  and  the  remote  object  at  such 
distances,  the  time  difference  is  on  the  order  of  10~16  seconds.  Because  the  coherence 
time,  tc,  of  a  conventional  laser  is  inversely  proportional  to  the  laser  bandwidth  [17], 
a  coherent  illumination  across  the  extent  of  the  object  is  easily  achieved  with  a  laser 
with  bandwidth  less  than  1GHz.  This  laser  bandwidth  is  easily  achieved  with  available 
technology  and  does  not  present  cost  or  design  concerns  for  a  proposed  system. 


Figure  3.3:  Geometry  for  Considering  Spatial  Coherence 

Because  polarization  sensitivity  is  to  be  added,  laser  illumination  characteris¬ 
tics  should  be  considered  during  system  design.  We  desire  statistically  independent 
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speckle  realizations  in  the  two  observed  channels.  Reflectivity  as  a  function  of  po¬ 
larization  may,  depending  on  the  object  surface,  provide  statistical  independence 
between  polarimetric  data  sets.  However,  further  system  design  should  be  adopted  to 
ensure  this  effect  is  more  readily  observed.  Simultaneous  illumination  with  two  laser 
beams  with  orthogonal  polarizations  and  phase  front  difference  (e.g.  tilt)  will  also  pro¬ 
duce  different  speckle  realizations  within  the  two  channels.  Certainly,  with  sufficient 
phase  differences  in  the  two  illuminating  beams,  different  realizations  of  reflected 
phase  are  created  at  the  object  surface.  Also,  at  low  light  levels,  photo-detection 
noise  will  dominate  the  detection  process  providing  statistically  independent  noise 
realizations  in  the  two  channels  because  of  separate  photo  detectors. 

3.2  Mathematical  Model 

The  illuminated  field  at  the  target  plane,  /( u),  is  spatially  coherent  across  the 
extent  of  the  target  (or  beam  extent).  The  reflected  field  is  observed  repeatedly  over 
many  laser  pulses  in  the  far-field  by  a  synthetic  array  of  optical,  “light-bucket”  de¬ 
tectors  without  the  aid  of  an  optical  lens.  K  statistically  independent,  laser  speckle 
patterns  are  transformed  through  post-processing  to  produce  noisy  object  autocorre¬ 
lations,  c4(x)- 


4(x)  =  |4-1{4(v)}|2,  (3.1) 

where  4(v)  is  the  fcth  frame  of  observed  laser  speckle  data  and  is  the  inverse 
Fourier  transform  performed  digitally  in  a  computer  via  a  Digital  Fourier  Transform 
(DFT).  The,  hth  frame  of  nonimaged  laser  speckle  data,  4(v),  is  modeled  by  a 
Fraunhofer  propagation  (JFA2)  with  a  mean  wavelength  A  and  propagation  distance 

A 


4(v)  =  |4A4/fc(u)}|-  •  A(v)  +  nfc(v),  (3.2) 
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where  A(v)  is  the  aperture  function  denoting  the  region  the  speckle  pattern  is  ob¬ 
served  and  n*;(v)  is  noise  encountered  in  the  detection  process  that  includes  pho¬ 
ton  noise,  read-out  noise,  and  noise  due  to  background  light.  Speckle  noise  is  in¬ 
herent  and  occurs  prior  to  detection  due  to  coherent  illumination  and  the  random 
phase  imparted  on  the  reflected  wavefront  by  the  rough  surface  of  the  target  (e.g. 
/(u)  =  |/(u)|  exp{j#(u)};  0(u)  is  modeled  as  uniformly  random  ~  U[— tt,  7t]  and 
independent).  Normally,  a  laser  line  filter  is  used  to  minimize  background  light. 

Recalling  Eqn.  12.21  the  mean  of  the  noisy  object  autocorrelations  is  related  to 
the  true  object  autocorrelation,  R0(x), 

E[  dk(x)  ]  =  b\h(x.)\2  +  c[i20(x)  *  |/i(x)|2],  (3.3) 

where  b  and  c  are  constants,  |/i(x)|2  is  the  known,  incoherent  point  spread  function 
(PSF)  of  the  detector  array,  i?G(x)  is  the  true  autocorrelation  of  the  object  intensity 
(o  =  |/|2),  and  *  denotes  convolution.  If  the  detector  array  is  uniform  and  no  zero- 
padding  is  used  in  the  DFT,  the  PSF  is  a  weighted  Dirac  delta  function,  <5(x),  and 
the  mean  of  the  transformed  data  simplifies  to 

E[  4(x)  ]  =  b8(x)  +  cR0(x).  (3.4) 

The  strength,  b  of  the  delta  function  in  EqifJ.i3.4iis  related  to  the  detector  array  and  the 
intensity  of  the  object  scene  (see  Appendix  [B]).  To  avoid  estimation  complexity  for  a 
single  pixel  of  the  data  image,  this  pixel  value  is  removed  prior  to  any  computation  and 
replaced  with  the  peak  of  the  estimated  object  autocorrelation  from  the  initial  guess 
or  old  object  estimate,  oold,  from  the  previous  iteration.  With  these  simplifications,  a 
measured  object  autocorrelation,  Ra(x)  is  obtained.  The  measured  autocorrelation, 
i2o(x),  is  defined  by: 
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(3.5) 


K 


R0{x)  =  K  1  5^4(x)[!  -  5(x)]  +  [°old(y)]2(5(x)- 


k= 1 


y 


The  measured  autocorrelation  is  the  average  of  the  observed  data,  4(x),  with  the 
central  image  pixel  modified  as  discussed  above. 

3.3  Polarimetric  Model 

The  polarimetric  data  in  this  research  is  obtained  via  a  two-channel  system: 
one  channel  is  polarization  insensitive  and  the  second  channel  is  observed  through 
polarizing  film  or  via  a  polarization  beam  splitter  as  detailed  in  Sec.  3.1.1.  Originally 
proposed  by  Strong  [43],  a  lumped  parameter,  p,  is  introduced  that  is  the  polarization 
ratio.  The  polarization  ratio  is  the  ratio  of  the  intensity  observed  in  the  polarized 
channel  to  the  intensity  observed  in  the  unpolarized  channel, 


(3.6) 


where  0  <  p  <  1.  The  polarization  ratio  is  essentially  the  projection  of  the  object 
intensity,  o,  as  viewed  through  the  polarizer,  or,  the  fraction  of  light  transmitted 
through  the  polarizer, 


0r(  y)  =  p(y)o(y). 


(3.7) 


For  this  development,  the  polarization  of  the  object  scene  (or  degree  of  polarization) 
is  not  estimated  nor  is  an  assumption  made  regarding  the  decomposition  of  the  re¬ 
flected  light  into  polarized  and  unpolarized  components  as  found  in  [28].  The  lumped 
parameter  p  is  not  useful  for  determining  the  scene’s  degree  of  polarization;  however, 


this  parameter  enables  us  to  relate  polarized  and  unpolarized  data  elements.  The 
polarization  ratio,  p  is  only  used  to  solve  for  the  unknown  object,  o. 
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Normally,  this  phase  retrieval  problem  can  be  characterized  by  one  equation 
and  two  unknowns  (Fourier  magnitude  and  phase).  By  adding  a  second  measure¬ 
ment  that  is  related  to  the  first,  the  problem  is  characterized  by  two  equations  and 
three  unknowns  (magnitude,  phase,  and  polarization  ratio).  The  second,  but  related 
measurement,  further  constrains  the  problem,  improving  search  performance.  Us¬ 
ing  Eqn.  13.4}  suppressing  the  energy  at  the  central  pixel,  b<5(x),  and  applying  the 
polarimetric  model,  the  Expected  Values  of  the  data  in  the  two  channels  are 


E[dk\*)\  =  [1-^(x)]E°(y)°(y  +  x)' 

y 

(3.8) 

E[d^\x)}  =  [l-5(x)]^op(y)op(y  +  x). 

(3.9) 

y 


Equations  13.81  and  3.9  are  used  in  Chapter  [TV]  to  develop  an  algorithm  to  solve  for 
the  unknown  object. 

3-4  Statistical  Model 

An  accurate  statistical  model  is  important  for  any  statistical-based  estimator, 
the  focus  of  this  research.  Previous  research  in  correlography  similar  to  the  application 
considered  here  has  not  explored  the  statistics  of  the  processed  data.  Much  of  the 
image-recovery  via  phase-retrieval  research  area  involves  gradient  search,  root-finding 
or  Fourier  transform-based  algorithms  where  a  “best-fit”  solution  is  found  through 
iterative  search  techniques  [12],  Often,  assumptions  about  the  remote  object  are 
used  to  make  the  computation  tractable  (e.g.  positivity  and  spatial  bound).  These 
approaches  do  not  employ  a  statistical  model  for  the  data  (e.g.  [11,24]).  Other 
research  efforts  provide  an  assumption  of  the  noise  statistics  for  the  processed  data 
(e.g.  Gaussian  [23,33]).  This  research  endeavors  to  better  characterize  the  noise 
statistics  of  the  processed  data.  First,  the  underlying  assumptions  of  the  object  are 
stated  and  assumed  similar  to  those  found  in  previous  research. 
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The  statistics  of  laser  speckle  due  to  polarized  thermal  light  are  well  known 
[16,17,19].  The  probability  distribution  function  (pdf)  of  instantaneous  intensity  for 
speckle  caused  by  object  phase  with  a  uniform  distribution  is  known  to  be  negative 
exponential.  The  intensity  captured  by  any  real  detector  will  be  an  integrated  in¬ 
tensity  for  a  finite  integration  time,  T,  and  measured  in  photocounts.  The  resulting 
density  function  is  the  negative  binomial  distribution  [17].  This  assumes  perfectly 
polarized  light.  Partially  polarized  or  unpolarized  light  will  have  a  different  result. 

However,  with  three  key  assumptions,  the  statistics  of  the  observed  laser  speckle 
intensity  may  be  analyzed: 

1.  The  amplitude  and  phase  of  the  reflected  field  are  statistically  independent. 
This  assumption  is  reasonable  since  the  amplitude  is  a  function  of  the  object 
reflectivity  and  the  phase  is  a  function  of  the  surface  roughness  or  height  profile. 
These  two  physical  elements  are  unrelated. 

2.  The  phase  is  spatially  independent  and  identically  distributed.  This  assumption 
is  essentially  that  the  sample  size  is  not  small  enough  to  produce  statistical 
dependence  due  to  nearness  of  the  sampled  points  being  related  by  similar 
roughness. 

3.  The  object  surface  roughness  is  modeled  as  a  uniform  random  variable  dis¬ 
tributed  between  (— vr,  tt)  and  statistically  independent.  This  assumption  is 
valid  for  most  surfaces  (other  than  mirrored  surfaces)  as  the  surface  is  consid¬ 
ered  “rough”  if  the  surface  height  profile  is  much  greater  than  the  wavelength 
of  the  illuminating  field. 

Even  if  the  distribution  of  the  observed  laser  speckle  intensity  deviates  from  the 
negative  binomial  due  to  partial  polarization  of  the  illuminating  beam,  a  key  result  is 
observed.  If  these  three  assumptions  are  valid,  the  expected  value  of  the  laser  speckle 
intensity  is  easily  shown  to  be  a  constant  (see  Appendix  [AT). 

Knowledge  of  the  average  value  of  the  laser  speckle  intensity  patterns  will  enable 
the  exploration  of  the  statistics  of  the  processed  or  transformed  laser  speckle  data. 


38 


Largely  due  to  this  expected  value  result,  Appendix  IQ  demonstrates  the  statistics  of 
the  transformed  data  is  well  approximated  by  the  exponential  distribution.  Because  of 
this  result,  an  ML  solution  was  investigated  for  solving  this  particular  phase  retrieval 
problem. 

It  will  be  assumed  that  each  frame  of  collected  laser  speckle  intensity  data  is 
statistically  independent.  This  is  accomplished  if  the  phase  at  the  object  surface  is 
statistically  independent  from  pulse  to  pulse.  This  is  reasonably  observed  if  minute 
geometry  changes  occur  from  pulse  to  pulse  that  produce  this  effect.  These  pulse- 
to-pulse  geometry  changes  may  include  target  jitter  due  to  movement,  target  surface 
deformation  due  to  compliant  structures,  laser  (source)  jitter,  etc.  Certainly,  this 
effect  is  a  system  design  consideration  as  each  frame  of  laser  speckle  must  be  statisti¬ 
cally  independent  for  the  result  detailed  in  Appendix  0  to  be  observed.  Additionally, 
it  is  assumed  the  simultaneously  observed  data  in  the  two  channels  (same  pulse  or 
frame)  as  collected  via  a  PBS  (S  and  P  channels)  are  statistically  independent  due 
to  surface  reflectivity  as  a  function  of  polarization  and/or  system  design  with  dual 
illumination.  Simultaneous  illumination  with  two  laser  beams  with  orthogonal  po¬ 
larizations  and  phase  front  difference  (e.g.  tilt)  will  also  produce  different  speckle 
realizations  within  the  two  channels.  Also,  at  low  light  levels,  photo-detection  noise 
will  dominate  the  detection  process  providing  statistically  independent  noise  realiza¬ 
tions  in  the  two  channels  because  of  separate  photo  detectors. 

It  will  be  further  assumed  the  processed  data  frames  are  also  statistically  inde¬ 
pendent.  It  is  assumed  the  mathematical  transformation  (magnitude  squared  of  the 
Fourier  transform)  operating  on  each  frame  does  not  create  statistical  dependence 
between  data  frames. 

With  exponential  statistics,  the  Signal-to- Noise  Ratio  (SNR)  is  easily  described. 
The  exponential  distribution  is  completely  described  by  the  first  moment  or  the  ex¬ 
pected  value.  Also,  the  expected  value  equals  the  standard  deviation.  For  this  reason, 
the  SNR  of  the  data  set  is  a  function  of  the  number  of  frames  collected,  K.  Because 
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the  K  noisy,  exponentially  distributed  autocorrelations  are  summed  together,  the 
final  distribution  of  the  average  autocorrelation  is  distributed  as  a  gamma  random 
variable  with  expected  value  equal  to  K/ R0  and  standard  deviation  equal  to  \[K / Ra. 
Therefore,  the  SNR  for  the  averaged  autocorrelation  computed  from  K  frames  of  data 
is  equal  to  \/~K. 

3.5  Maximum  Likelihood  Approach 

With  laser  speckle  intensity  data  transformed  to  noisy  autocorrelations  and  the 
noise  modeled  with  an  exponential  probability  density  function,  an  ML  solution  is 
investigated.  Only  a  single  channel  of  unpolarized  data  is  considered. 

The  observed  data  is  the  autocorrelation  of  the  unknown  object  corrupted  by 
exponential  noise.  Many  statistically  independent  realizations  are  observed  and  col¬ 
lected.  It  is  assumed  the  noise  at  each  sampled  point  in  the  data  image  is  statistically 
independent  from  all  other  points.  The  joint  probability  density  function  (pdf)  for  a 
single  frame  of  data  is 


PD(d(x))  =  (3-10) 

X  °  '  ' 

The  expected  value,  Ra,  is  defined  by 

^o(x)  =  J^o(y)o(y  +  x),  (3.11) 

t 

where  o(t)  is  the  unknown  object  intensity  and  y  and  x  are  two  dimensional  co¬ 
ordinate  vectors.  We  will  also  assume  each  frame  of  observed  data  is  statistically 
independent  from  all  other  collected  data  frames.  For  K  frames  of  data,  the  joint  pdf 

is 


PD{di{x), ...,  dK{x))  =  JJU 

k  x 


1 

R0(x)C 


4fc(x) 

floO)  . 


(3.12) 
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Forming  the  log-likelihood  function,  L,  by  taking  the  natural  logarithm  of  the  joint 
pdf  yields 


L  =  ~KY,lnR  (3-13) 

Taking  the  partial  derivative  of  L  with  respect  to  the  object  produces  the  next  equa¬ 
tion. 


dL  _  k'S~'  ^  9R0(x)  4(x)  dR0(x) 

do( y)  “  '  4“  ^x)  ^  k  x  ^o(x)  do(y)  1  '  ■ 

Next,  the  partial  derivative  of  i?c(x)  with  respect  to  the  object,  o(y),  is  computed. 


<9R0(x) 
do(  y) 


o(y  +  x)  +  o(y  -  x) 


Substituting  this  result  into  dL/do  yields 


(3.15) 


=  ~K  Y  Ro  X(x)  [o(y  +  u)  +  o(y  -  u)] 

+  YY1  R° 2(x)4(x)  [o(y  +  u)  +  o(y  -  u)] .  (3.16) 

k  x 

Setting  this  equal  to  zero  and  solving  for  o  maximizes  the  function  L  with  respect  to 
the  object,  o.  Since  a  closed  form  solution  for  o  is  not  feasible,  an  approach  similar 
to  the  Richardson-Lucy  (RL)  algorithm  used  in  deconvolution  [29,  34]  is  employed. 
The  ratio  of  the  negative  to  positive  parts  of  this  function  is  formed  and  o  is  solved 
iteratively  from  an  initial  guess.  First,  define  the  average  value  of  the  observed  data, 
.D(x),  as 


D(x)  =  ^^4(x)- 

k 


(3.17) 
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The  RL  approach  to  iteratively  solve  for  o  is  then  given  to  be: 


onew(  y)  =  o°ld(  y) 


£M 

RiM 


RoM 


-k  o(x) 
-k  o(x ) 


+ 

+ 


DM 

RoM 


RoM 


*  o(x) 

*  o(x)  ’ 


(3.18) 


where  *  is  correlation  and  *  is  convolution.  This  algorithm  is  problematic  in  two 
respects:  (1)  The  algorithm  does  not  naturally  constrain  the  object  magnitude,  and 
allows  for  o  to  grow  without  bound  each  iteration,  and  (2)  the  algorithm  presents 
numeric  challenges  due  to  division  by  the  squared  term  Ra  which  is  the  current 
estimate  of  the  object  autocorrelation.  This  ML  approach  utilizes  a  more  correct 
statistical  model  than  previously  published  algorithms;  however,  it  does  not  provide 
a  useable  algorithm  as  evidenced  by  analysis  with  computer  simulation.  The  general¬ 
ized  expectation-maximization  technique  provides  a  broader  and  more  powerful  ML 
estimate  compared  to  the  ML  technique  described  above. 


3.6  Expectation- Maximization  Approach 

The  EM  algorithm  was  systematically  defined  and  convergence  proved  in  the 
seminal  paper  by  Dempster,  Laird  and  Rubin  (DLR)  [9].  The  DLR  paper  coalesced 
previous  research  and  journal  papers  treating  this  generalized  approach  to  developing 
ML  estimates.  The  EM  technique  has  wide  applicability  with  desirable  convergence 
properties.  This  approach  provides  a  powerful  tool  for  solving  problems  involving 
missing  or  incomplete  data  or  direct  access  to  the  data  necessary  to  estimate  the 
required  parameters  is  impossible  [30].  Although  convergence  properties  have  been 
revisited  since  the  original  paper  [45] ,  the  EM  algorithm  does  guarantee  convergence 
to  a  local  maximum  as  the  likelihood  function  is  increased  at  every  iteration.  The 
EM  algorithm  is  guaranteed  to  be  stable  and  converges  to  an  ML  estimate  [30].  The 
EM  approach  is  well-suited  for  the  phase  retrieval  problem  presented  in  this  research. 
The  EM  approach  has  been  used  in  a  related  phase  retrieval  problem  [39]. 

We  wish  to  estimate  an  object,  o(x),  from  many  statistically  independent  real¬ 
izations  of  observed  data  produced  by  a  random  process.  However,  the  observed  data 
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are  corrupted  by  noise.  From  observed  data,  noisy  autocorrelation  data  are  produced 
via  a  post-processing  transformation  operation.  Each  point  in  the  resulting  autocor¬ 
relation  data  is  an  identically  distributed  and  statistically  independent  exponential 
random  variable.  This  processed  data  provides  an  incomplete  view  of  the  parame¬ 
ter  to  be  estimated  and  is  called  incomplete  data.  The  unobserved  data  containing 
the  required  information  is  called  the  complete  data,  4( y,  x).  The  complete  data  is 
related  to  the  incomplete  data,  4(x)  by 


4(x) 


2 


^2  <4(y.  x)eJ"(y) 


(3.19) 


where  the  complete  data  is  multiplied  by  a  uniformly  distributed  phasor,  9(y),  summed 
over  all  values  of  y,  and  then  a  magnitude  squared  operation  is  performed.  The 
subscript  k  denotes  the  data  frame  and  x  and  y  are  each  two-dimensional  spatial 
variables.  This  operation:  random  phasor  sum,  magnitude  squared  operation  results 
in  a  quantity  with  an  exponential  distribution.  The  incomplete  data  is  known  to 
be  exponentially  distributed  with  mean  equal  to  the  autocorrelation  of  the  desired 
object.  The  random  phasor  sum,  magnitude  squared  seems  to  be  a  natural  choice  of 
complete  data  that  leads  to  exponentially  distributed  incomplete  data. 


£[4(x)]  =  Y  °(y)°(y  +  x)  (3.20) 

y 

Because  the  the  mean  of  the  incomplete  data  is  known,  the  complete  data  is  chosen 
to  have  the  following  property: 


E[dl(  y> x)]  =  o(y)o(y  +  x).  (3.21) 

The  complete  data  is  selected  in  this  manner  due  to  the  following: 
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£[4(x)]  =  E 


^2^dk(y:it)dk(y i  +x) exp{j[0(y)  -  0(yi)]} 
y  yi 


=  ££  JB[4(y,x)4(yi  +x)]JE[exp{j[0(y)  -  ^(yx)]}] 
y  yi 


=  ^2^2E[dk(y,  x)4(yx  +x)]5(y-y1) 
y  yi 


=  J2E[dl(y,^)] 

y 


=  ^o{y)o{y  +  *) 
y 


(3.22) 


Because  the  phase  term,  9(y)  is  uniformly  distributed  [— tt,  7t]  and  independent,  the 
complete  data  can  be  of  any  distribution  or  non-random  and  yield  the  correct  mean 
for  the  incomplete  data.  The  distribution  of  the  complete  data  can  then  be  chosen  in 
the  most  advantageous  manner.  The  square  of  the  original  complete  data  is  chosen 
to  be  the  new  complete  data  of  interest  and  a  Poisson  random  variable. 


4(y,x)  =d\{  y,x) 


(3.23) 


P 


4(y,  x),  ...,dK(y,  x) 


nnn 


[o(y)o(y  +  x)]dfc(xx) 
4(y>x)! 


exp  [  -  o(y)o(y  +  x)] 

(3.24) 


The  log-likelihood  function  of  the  complete  data,  Lqd ,  is  found  by  taking  the  natural 
log  of  the  probability  mass  function  in  Eqn.  13.241 
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(3.25) 


Lcd(o)  =  EEE  4(y,  x)  log[o(y)o(y  +  x)]  -  [o(y)o(y  +  x)] 

k  x  y  ^ 

-  log[4(y,x)!]| 

3.6.1  Expectation  Step.  The  expectation  step  of  the  EM  algorithm,  or 
Q- Function  (Q),  is  defined  as  the  expectation  of  the  complete  data  log- likelihood 
function  conditioned  on  the  old  estimate  of  the  object  from  the  previous  iteration, 
oold(y),  and  the  incomplete  data,  4(x), 

Q(o|o°ld,4(x))  =  E[Lcd(o)  |  oold,4(x)].  (3.26) 


Q(o\o°ld)  =  ^'52'52E[dk(y,x)\o°u,dk(x)]  -  log  [o(y)o(y  +  x)] 

k  x  y 

-  K  [°(y)°(y  +  x)]  -  a.t.  (3.27) 

x  y 

where  K  is  the  total  number  of  frames  and  A.T.  denotes  another  term  not  a  function 
of  the  object,  o.  Next,  the  conditional  expectation  of  the  complete  data  given  the 
incomplete  data,  /i,  is  computed: 


/i(o°ld,  4(x))  =  E  [4(y,  x)|oold,  4(x)] .  (3.28) 

This  is  often  the  most  difficult  step  in  computing  the  EM  algorithm,  ft  may  be 
difficult  to  solve  for  the  mean  of  the  conditional  density  function.  This  was  attempted 
for  Eqn.  3.281  Using  Bayes’  rule,  we  define  the  conditional  probability  mass  function 
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(pmf)  of  the  complete  data  given  the  incomplete  data  (for  a  specific  frame  k,  and 
specific  spatial  variables  y  and  x)  as 


P[dk( y,x)  j  4(x)] 


p[4(x)  1  4(y,x)]P[4(y,x)] 
p[4(x)] 


(3.29) 


The  probability  mass  function  of  the  complete  data  is  specified  and  the  prob¬ 
ability  density  function  of  the  incomplete  data  is  known.  The  probability  density 
function  of  the  incomplete  data  given  the  complete  data  can  be  assumed  to  be  expo¬ 
nentially  distributed;  therefore  can  be  specified  by  only  its  mean.  The  mean  of  this 
conditional  density  is  found  by 


£[4(x)  |  4(y0,x)]  =  EE  E[4(y,x)4(yi  +  x)  |  4(yO)x)]£[exp{j[0(y)  -  0^)]}] 

y  ?/i 

=  EE  E[4(y,x)4(yi +x)  |  4(y0>x)]<J(y-y1) 

y  2/1 

=  5^-E[4(y,x)  |  4(y0, x)] 

y 

=  Y  °(y)°(y  +  x)  -  o(y0)o(y0  +  x)  +  4(y0,  x).  (3.30) 

y 

From  this  result  the  probability  density  function,  pj),  of  the  incomplete  data  condi¬ 
tioned  on  the  complete  data  is  given  as 


Pz?[4(x)  |  4(y,x)] 


exp 


_ -<4(x) _ 

E y  °(y)°(y+x) -o(y0 )o(y0+x)+dfc (y0 ,x) 


E„°(y)o(y  +  x)  -  o(y0)o(y0  +x)  +  4(y0,x) 


(3.31) 


Returning  to  Bayes7  rule,  this  specifies  the  probability  mass  function  of  the  incomplete 
data  given  complete  data. 
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P[4(y,x)  |  4(x)]  = 


E„°(y)o(y  +  x) 


lE, °(y)°(y  +  x)  -  °(y0)°(yo  +  x)  +  4(y0, x) 


■  exp 


exp 


— 4(x) 


lE y  o(y)o(y  +  x)  -  o(y0)o(y0  +  x)  +  dk(y0,  x) 


4fx) 


Ey  o(y)o(y  +  x) 

exp[-o(yo)o(y0  +  x)] 
4(y0,x)! 


o(y0)o(y0  +  x) 


<4(  y0.x) 


(3.32) 


This  pmf  is  difficult  to  integrate  in  order  to  determine  the  mean.  However, 
approximations  were  attempted  in  order  to  find  functions  easily  integrable.  With 
some  reasonable  approximations,  this  pmf  may  be  simplified  into  a  recognizable  form 
where  the  mean  may  be  surmised.  However,  these  approximations  led  to  an  algorithm 
solution  that  did  not  properly  converge.  Adding  a  second,  but  related  polarimetric 
channel,  only  makes  this  more  difficult.  Without  a  good  solution  for  the  conditional 
expectation  of  the  complete  data,  /i(oold,  4(x)),  the  EM  approach  using  the  exponen¬ 
tial  noise  model  was  not  pursued  further. 


3.7  Poisson  Statistics 

The  exponential  model  was  unsuccessfully  explored  as  detailed  above.  There¬ 
fore,  a  different  statistical  model  was  explored.  It  is  known  the  measured  autocor¬ 
relation  is  not  corrupted  by  Poisson  noise;  however,  successful  algorithms  have  been 
developed  based  on  the  Poisson  model  for  this  specific  problem  [38,39].  Algorithms 
based  on  the  Poisson  model  for  this  type  problem  are  reasoned  to  minimize  the  I- 
divergence  measure  for  data  with  signal-dependent  noise  [38].  Both  the  Schulz  and 
Snyder  [38]  algorithm  and  the  Schulz  and  Voelz  algorithm  [39]  demonstrate  attrac¬ 
tive  properties  and  meaningful  results.  Also,  the  Poisson  distribution  is  characterized 
only  by  its  mean,  similar  to  the  exponential  distribution.  The  Poisson  model  is  again 
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chosen  for  this  research.  The  ML  phase  retrieval  algorithm  using  Poisson  statistics  de¬ 
scribed  by  Schulz  and  Snyder  [38]  is  re-derived  in  Appendix  JJl  using  an  EM  approach. 
The  resulting  iterative  estimator  for  the  unknown  object  is 


°new(y0)  = 


2S" 


Oo*  *  +  0°U  *  R' 


R°}d 


R°ld 


(y. 


0  )i 


(3.33) 


where  *  is  correlation,  *  is  convolution,  onew  is  the  new  object  estimate,  oold  is  the  old 
object  estimate  from  the  previous  iteration,  R0  is  the  measured  autocorrelation,  i?°ld 
is  the  autocorrelation  formed  from  the  old  estimate  of  the  object  from  the  previous 
iteration,  and  S “ew  is  the  two-dimensional  sum  of  the  new  object  estimate  computed 
by  Eqn.  ID. 15  in  Appendix  [Dl 

This  iterative  algorithm  derived  via  EM  technique  should  be  similar  or  exactly 
equal  to  the  Schulz  and  Snyder  ML  algorithm  [38]  (see  Eqn.  |2.3) .  By  inspection, 
the  EM  derived  algorithm  only  differs  from  the  Schulz  and  Snyder  ML  algorithm  by 
a  scale  factor  (see  Appendix  |D).  This  result  is  provided  here  for  comparison  to  a 
multi-channel  approach  detailed  in  Chapter  HVl 
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IV.  Polarimetric  Algorithms 

This  chapter  describes  two  EM  algorithms  for  use  with  polarimetric  data.  A  maximum 
likelihood  estimator  is  formed  to  process  polarimetric  data  related  to  the  unknown 
object’s  autocorrelation.  Two  different  system  models  are  explored:  (1)  a  two-channel 
system  and  (2)  a  dual-channel  system.  The  difference  in  the  system  models  is  the 
polarization  configuration  for  the  two  data  channels.  The  EM  technique  employed 
here  follows  closely  with  the  clear  development  of  the  generalized  EM  algorithm  found 
in  [37,39], 

The  EM  algorithm  technique  allows  for  generalized  developments  such  as  an 
extension  to  maximum  a  posteriori  (MAP)  estimators  as  shown  by  Dempster,  Laird 
and  Rubin  [9].  In  this  chapter,  the  EM  method  is  extended  to  a  MAP  approach  for 
the  two-channel  system  with  the  introduction  of  a  prior  distribution  on  the  polariza¬ 
tion  parameter,  p.  Lastly,  a  statistical-based  stopping  criteria  is  provided  for  timely 
stopping  of  the  iterative  algorithms. 

4-1  Two-Channel  Algorithm 

A  two-channel  system  observes  unpolarized  and  polarized  data  in  two  channels. 
Channel  one  produces  the  unpolarized  data  set  and  channel  two  produces  the  polar¬ 
ized  data  set.  As  presented  in  Chapter  HU  there  are  multiple  approaches  to  system 
design  to  produce  the  polarized  and  unpolarized  data  sets.  However,  this  correlog- 
raphy  method  requires  statistical  independence  in  the  two  channels,  an  important 
design  consideration.  The  following  development  delineates  an  EM  algorithm. 

4-1.1  Incomplete  Data  Model.  The  incomplete  data,  c4(x),  is  a  series  of 
statistically  independent,  noisy  autocorrelations  of  the  original  object  transformed 
from  laser  speckle  observations.  With  an  assumption  about  the  target  scene’s  surface 
roughness,  the  observed  laser  speckle  data  is  known  to  be  statistically  distributed 
as  Negative  Binomial  due  to  the  doubly  stochastic  process  of  speckle  generation  and 
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photoelectric  detection  [17].  However,  the  transformation  of  the  laser  speckle  data 
into  noisy  autocorrelations  (see  Eqn.  13. lj)  produces  entirely  different  statistics  and  can 
be  approximated  by  an  exponential  distribution  with  mean  equal  to  the  true  object 
autocorrelation  (see  Appendix  A).  The  exponential  model  was  unsuccessfully  explored 
as  detailed  in  Chapter  II1I1  The  Poisson  model  is  chosen  for  the  development  of  a 
multi-channel  polarimetric  algorithm.  This  development  is  a  polarimetric  extension 
to  the  Schulz  and  Snyder  algorithm  for  recovering  images  from  correlation  data  [38] 
developed  with  the  EM  technique.  The  single  channel,  non-polarimetric  variant  is 
developed  for  completeness  in  Appendix  [EQ 

4-1-2  Complete  Data  Model.  The  complete  data  is  postulated  to  be  statisti¬ 
cally  independent  variates,  c4( y,  x),  distributed  as  Poisson,  related  to  the  incomplete 
data  and  with  expected  values: 


41) (x)  =  E5fc1)(y>x)«  t4-1) 

y 

E$k\  y,x)]  =o(y)o(y  +  x),  (4.2) 

42) (x)  =  E42)(y,x),  (4-3) 

y 

E$k\  y>x)]  =  oP(y)°p(y  +  x),  (4-4) 


where  E[-\  is  the  expected  value  operator,  x  and  y  are  two-dimensional  spatial  vari¬ 
ables,  o  is  the  unknown  object,  op  is  the  unknown  polarized  object  as  viewed  through 
the  polarization  analyzer,  and  (1)  and  (2)  indicate  channel  number.  Assuming  sta¬ 
tistical  independence  between  the  two  channels  and  each  frame,  the  complete  data 
log-likelihood,  Lqd,  is 
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K  K  , 

lCd(o)  =  1 41}(y>x)  log[°(y)°(y +  x)]  +  ^2)(y>x)  logMyK(y  +  x)] 

k= 1  1=1  x  y 

-  o(y)o(y  +  x)  -  Op(y)op(y  +  x)  +  O.T.j,  (4.5) 

where  k  and  /  are  independently  indexed  frame  numbers  and  O.T.  represents  other 
terms  not  a  function  of  o  or  p.  These  terms  are  eliminated  in  the  maximization  step; 
therefore,  they  are  ignored. 

With  this  formulation  of  the  log-likelihood  function,  the  EM  solution  for  the 
unknown  object  o  degenerates  into  a  solution  as  function  of  the  data  observed  only 
in  the  unpolarized  channel  (e.g.  onew  is  not  a  function  of  pnew).  This  degeneration 
of  the  two-channel  polarimetric  system  also  occurs  in  the  imaging  case  as  detailed 
by  Strong  [43].  To  overcome  this  degenerative  solution  and  use  all  of  the  data  from 
both  channels,  Strong  proposed  a  departure  from  the  EM  technique  by  using  an 
old  estimate  of  the  polarization  ratio,  pold,  for  calculating  onew.  Strong  successfully 
used  this  substitution  with  reasonable  results  but  produced  an  algorithm  no  longer 
characterized  as  EM.  This  research  explored  a  similar  approach  with  satisfactory 
results  for  o;  however,  the  estimated  values  for  p  become  non-physical  and  are  often 
estimated  to  be  much  larger  than  one.  To  overcome  this  difficulty  and  maintain 
the  EM  algorithm  technique,  it  is  proposed  to  include  a  prior  distribution  on  the 
polarization  ratio,  p,  similar  to  the  imaging  case  found  in  [26].  The  introduction  of 
the  prior  distribution  extends  the  EM  algorithm  to  a  MAP  estimator  vice  an  ML 
estimator. 

The  polarization  ratio  is  known  to  be  positive  but  also  less  than  or  equal  to 
one;  however,  a  uniform  density  function  (e.g.  p  ~  f/[0, 1])  is  not  helpful.  It  can  be 
reasoned  that  the  polarization  ratio  of  a  random  scene  is  less  likely,  on  average,  to 
produce  a  p  equal  to  zero  or  one  and  more  likely,  on  average,  to  produce  a  p  near  0.5. 
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Therefore,  a  normal  distribution,  with  density  function  f(p),  centered  at  0.5,  is  an 
ideal  prior, 


(4.6) 


where  s  is  a  scaling  parameter,  a  is  an  arbitrary  shape  parameter,  and  n  is  a  posi¬ 
tive,  even  integer.  However,  this  nth  order  distribution  adds  unwanted  computational 
complexity  to  the  estimator.  Consequently,  the  selection  of  a  meaningful  prior  dis¬ 
tribution  must  be  as  simple  as  possible  to  keep  the  problem  tractable.  Additionally, 
any  algorithm  development  must  enforce  positivity  for  both  o  and  p.  Therefore,  the 
simple  exponential  distribution  with  a  mean,  A  =  1/2,  is  chosen.  This  choice  of  A  is 
arbitrary  but  found  to  perform  reasonably  well.  The  exponential  distribution  with 
density  function,  f(p),  constrains  p  >  1  as  less  likely  compared  to  0  <  p  <  1, 


(4.7) 


This  distribution,  while  not  entirely  informative,  does  enforce  positivity,  and 


constrains  large  values  of  p  to  be  less  likely.  This  simple  selection  provides  these  two 


properties  and  enables  analytic  solutions  for  o  and  p  estimates  using  both  data  chan¬ 
nels.  With  this  selection  of  a  prior  constraining  p,  the  complete  data  log-likelihood 
function  becomes, 


K  K 


LCd(o,p)  =  \  41)(A’x)  log[°(y)°(y +  x)]  +  c?2)(y>x)  log[op(y)op(y  +  x)] 


k=  1  1=1  x  y 


k=  1  1=1 


52 


4-1.3  Expectation  Step.  For  the  Expectation  Step  (E-Step),  the  conditional 
expectation  of  the  complete  log-likelihood  is  taken: 


Q(o,p)  =E[Lcd\  dk(x),o  =  oold,p  =  p°,d ]  =  Eold[LCD\ 


IS  IS 


k=l  1=1  x  y 


+  £°"[tij'2)(y,x)]  log[op(y)Op(y  +  x)] 


(4.9) 


The  conditional  expectation,  Eold,  is  conditioned  on  the  incomplete  data  and  old 
estimates  of  the  object,  o,  and  parameter,  p.  The  conditional  expectation  of  the 
complete  data  is  often  the  most  difficult  step  of  the  EM  process.  Fortunately,  by 
choosing  the  Poisson  model,  the  form  of  the  conditional  expectation  is  provided  by 
Shepp  and  Vardi  [42], 


where  the  autocorrelations  formed  by  the  old  estimate  of  the  objects  o  and  op  are 
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(4.12) 


K“(x)  =  5>°ld(y)°°ld(y  +  *)• 

y 

r°poW  =  <d(y)°;ld(y  +  x)-  (4-13) 

y 

4-1-4  Maximization  Step.  The  Maximization  Step  (M-Step)  is  performed 
by  maximizing  Q(o,p )  of  Eqn.14.91  for  the  unknown  variates,  o  and  p.  The  Q-function 
is  maximized  by  finding  the  zero  of  the  Erst  partial  derivatives. 


9Q  =  V  V  V  /  i  Myp  -  x>x)  ,  M2(y0,x)  /i2(y0  -  x,x) 

do(y0)  hlkr\  °(yo)  °(yo-x)  o(y0)  o(y0) 

-  °(y0  +  x)  -  °(y0  -  x)  -  p(y0)°P(yo  +  x)  -  p(yoK(y0  -  x) 


=  o. 


(4.14) 


Solving  for  o(y0)  yields  the  new  estimate  for  o, 


°new(y0)  = 


1 


2 k2[s™  +  ^rrew( y0)]  t: 

+  h2(y0,x)  +  h2(y0  -x,x) 


K  K  , 

EEE  hi(y0ix)  +  My0  -x,x) 

k= 1  Z=1  x  ^ 


(4.15) 


where 


S0new  =  5>n6W(x)’ 

X 

5r=E°r(x)- 

X 


(4.16) 

(4.17) 
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Equation  14.151  can  be  simplified  by  evaluating  the  summation  terms: 


new/  X  =  ^l(y0)  +^2(y0) 

1  oJ  2[5'“ew  +  5“wpnew(y0)]  ’ 


(4.18) 


where 


^i(y0)  =  °old 
d'2(y0)  =  <d 


1  1 

o 

0^ 

K 

* 

1 _ 1 

(y0) 

+  oold 

oold*  Ro 

old  Rpo 

UP  E>old 

1  Lpo 

(y0) 

+  oold 

i 

old  RpO 

°P  *  7?old 

1  Lpo 

(yo). 

(yo). 


(4.19) 

(4.20) 


*  denotes  correlation,  and  *  denotes  convolution.  Note,  if  p  =  0  everywhere,  the 
second  channel  vanishes  and  the  above  solution  is  equivalent  to  the  single  channel 
solution  [38]  solved  via  ML  technique.  However,  if  p  >  0  for  any  point  in  the  scene, 
the  new  terms  assert  themselves  to  yield  the  proposed  multi-channel  algorithm. 

The  solution  for  onew  is  a  function  of  the  polarization  ratio,  pnew,  also  to  be 
estimated.  Therefore,  the  Q-function  will  be  maximized  for  the  parameter,  p.  This 
is  accomplished  by  also  finding  the  zero  of  the  first  partial  derivative, 


dQ 

My0) 


My0,x) 
P(  Yo) 


h2(y0-x,x) 

p(y0) 


o(y0K(yo  +  x)  -  o(y0)op(y0  -  x)  -  2 


=  0 


(4.21) 


Solving  for  p(y0)  yields  the  new  estimate  for  p, 


Pnew(Yo) 


^(y0) 

2[onew(y0)S'“ow  +  N }  ’ 


(4.22) 
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where  N  is  the  total  number  of  pixels  in  the  two-dimensional  spatial  vector,  x.  Plug¬ 
ging  Eqn.  14.221  back  into  Eq.  (4. 181  yields  a  quadratic  equation  from  which  the  positive 
root  is  selected.  This  provides  the  final  equation  for  estimating  onew, 


'(y0)  = 


Sp7^i  -  2 NS"™  +  [(2N  S°ew  -  S^Ti)2  +  8i\hS“ewS“w(^i  +  ^2)}1/2 


po 


po 


4S1" 


v  Onew 
upo 


(4.23) 


The  estimate  for  onew  is  a  function  of  old  estimates  and  the  data  in  both  polarized 
and  unpolarized  channels.  To  form  this  estimate,  both  the  scaling  constants  S'“ew  and 
Splw  must  be  computed.  First,  Q(o,p)  is  maximized  for  op  in  order  to  find  the  sum 
of  the  new  polarized  object  estimate,  Sp™, 


dQ 

dop(  y0) 


J/^2(yO)X)  .  At2(yo  — x;x)  /  ,  \  /  a 

— °>(yo+x>-^o-x>} 


k=  1  1=1 


=  0. 


(4.24) 


Solving  this  equation  for  S'"®"  and  summing  both  sides  of  the  equation  yields  a  solution 
for  the  scaling  factor  as  a  function  of  old  estimates  and  the  data  in  the  polarized 
channel, 


<?n 


'po 


(4.25) 


Finally,  the  sum  of  the  new  object  estimate,  S'®ew,  is  solved  for  by  summing  both  sides 
of  Eqn.  14.181, 


Qi  new  _ 


-5^[*i(y0)  +  *2(y0)]  -  (S', 


new  \  2 
po  ) 


y  o 


1/2 


(4.26) 
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4-2  Dual-Channel  Algorithm 

In  this  section,  a  dual-channel,  EM  algorithm  is  considered.  A  dual-channel 
system  observes  polarized  data  in  two  orthogonal  channels.  Employing  a  PBS  in  a 
non-imaging  LADAR  system  provides  two  channels  of  orthogonally  polarized  data. 
A  PBS  may  provide  polarized  data  with  minimal  light  loss  as  compared  to  traditional 
polarizers.  Also,  this  correlography  method  requires  statistical  independence  in  the 
two  data  channels  to  be  completely  effective;  a  requirement  aided  by  engineering  a 
system  with  a  PBS. 

4-2.1  Incomplete  Data  Model.  The  incomplete  data  is  the  observed  data 
with  the  object  obscured  by  noise  and  detection  limitations.  The  incomplete  data, 
d[^(x)  and  d[2\x),  is  a  set  of  statistically  independent,  noisy  autocorrelations  of 
the  polarized  object  transformed  from  laser  speckle  observations.  The  Poisson  noise 
model  is  again  chosen;  identical  to  the  two-channel  algorithm  with  the  same  ratio¬ 
nale.  The  expected  value  of  the  incomplete  data  is  modeled  as  autocorrelation  of  the 
polarized  object. 


Eldk\x)}  =  ^opl(y)opl(y  +  ^)  (4.27) 

y 

^K(2)(x)]  =  5^0p2(y)op2(y  +  X)  (4.28) 

y 

4-2.2  Complete  Data  Model.  The  complete  data  is  postulated  to  be  statisti¬ 
cally  independent  variates,  dk( y,  x),  distributed  as  Poisson,  related  to  the  incomplete 
data  and  with  expected  values: 
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(4.29) 


41}W  =  J341}(y.x)» 

y 

£[41}(y>x)]  =  Opi(y)opi(y +  x),  (4.30) 


42)(x)  =  5^42)(y’x)’  (4-31) 

y 

^[42)(y>x)]  =  Op2(y)oP2(y +  x).  (4.32) 

Assuming  statistical  independence  between  the  two  channels  and  each  frame,  the 
complete  data  log-likelihood,  Lcd ,  is 


K  K 

Lcd(o )  =  S  5Z{41}(y> x)  log[opi(y)opi(y  +  x)]  +  dj2)( y,  x)  log[op2(y)op2(y  +  x)] 

k= 1  Z=1  a:  j/ 

-  Opi(y)opi(y  +  x)  -  op2(y)op2(y  +  x)  +  O.T.},  (4.33) 

where  O.T.  represents  other  terms  not  a  function  of  o  or  p.  These  terms  are  eliminated 
in  the  M-Step;  therefore,  are  ignored. 

In  the  two-channel  case,  a  prior  distribution  on  p  was  introduced  to  avoid  a 
degenerate  solution.  However,  with  the  dual-channel  case,  a  prior  distribution  is 
not  required.  A  closed-form  solution  for  o  and  p  is  calculated  without  deviating 
from  the  EM  methodology.  Considering  orthogonality  of  the  two  channel  data,  the 
development  simplifies  the  number  of  unknowns.  Because  the  polarization  response  in 
channel  one  is  orthogonal  to  the  polarization  response  in  channel  two,  the  polarization 
ratios  of  the  two  channels  are  related  by: 
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My)  = 1  -pi(y)- 


(4.34) 


This  permits  the  dropping  of  the  channel  subscript  for  the  polarization  ratio  in  sub¬ 
sequent  equations  and  simplifies  the  polarized  object  model: 


Opi(y)  =p(y)o(y),  (4.35) 

Op  2(y)  =  [1  -  p(y)]o(y).  (4.36) 

4-2.3  Expectation  Step.  For  the  Expectation  Step,  the  conditional  expecta¬ 
tion  of  the  complete  log-likelihood  is  taken: 


Q(o,p)  =  E[Lcd\  4(x),o  =  oold, p  =  pold]  =  E°ld \LCd} 

K  K  , 

=  EEEE  ^[^(y,*)]  log[°pi(y)°pi(y  +  x)] 

k=l  1= 1  x  y  '• 

+  Eold[d}2)(y,x)]  log[op2(y)op2(y  +  x)] 

-  Opi(y)opi(y  +  x) +  op2(y)op2(y  +  x)|.  (4.37) 

The  conditional  expectation,  Eold,  is  conditioned  on  the  incomplete  data  and  old 
estimates  of  the  object,  o,  and  parameter,  p.  By  choosing  the  Poisson  model,  the 
form  of  the  conditional  expectation  is  known, 


My 


0)  - 


=  E° 


[41)(yo.x)]  =  o-(y0)0;?(y0  +  x) 


djpx) 

Kpi(x) 


(4.38) 
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My0,x)  =  £old$2)(y0,x)]  =  <2d(yo)<2d(yo  +  x) 


42)(x) 

*&(x)’ 


(4.39) 


where  the  autocorrelations  formed  by  the  old  estimate  of  the  polarized  objects,  op\ 
and  op 2,  are 


£l(x) 

=E°;'?(yH“(y+x), 

y 

(4.40) 

'op2(x) 

=  Y  o;2d(y)<2d(y  +  x). 

(4.41) 

4.2.4  Maximization  Step.  Again,  the  Maximization  Step  is  performed  by 
maximizing  Q(o,p)  of  Eqn.  14.371  for  the  unknown  variates,  o  and  p.  The  Q-function 
is  maximized  by  finding  the  zero  of  the  first  partial  derivatives. 

dQ  =  v  V  V  /  ^(yQ’x)  1  hi(y0  -  x,  x)  /i2(y0,x)  /i2(y0  -  x,  x) 

<9o(y0)  °^)  o(y0) 

-  2p(y0)p(y0  +  x)o(y0  +  x)  -  2[1  -  p(y0)]  [1  -  p(y0  +  x)]o(y0  +  x) 

=  0.  (4.42) 

Solving  for  o(y0)  yields  the  new  estimate  for  o, 


K  K  , 

/^i(yo.x)  +  /ii(y0  -x,  x) 

k  =  l  1=1  X  ^ 

(4.43) 


°new(y0)  = 


2 k2  [pnew(y0)^;i  +  [1  -  pnew(y0)]^;i]  “  jz 

+  /^2(y0;x) +  /i2(y0 -x,x)  l, 
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where 


=  5>rr(y)> 

X 

(4.44) 

s™2  =  5>;r(y)- 

(4.45) 

X 


Equation  14.431  can  be  simplified  by  evaluating  the  summation  terms: 


new/  }  =  _ dq(y0)  +  T2(y0) _ 

°  2  [p™ (y0)S%[  +  [1  -  Pnew(y0)]^;2]  ’ 


(4.46) 


where 


*i(y0)  =  <d 

^2(y0)  =  <2 


ou  Rppi 

>1  *  Dold 
■n'opl 

old  ,  -^op2 

°p2  Dold 
Jxop2 


(y0)  +  <i 
(y0)  +  <2 


oid  jypj 

°pl  Dold 

^opl 

old  .  Rop2 
0P2  Dold 

Jxop2 


(y0)> 

(y0)- 


(4.47) 

(4.48) 


.Ropi  and  -Rop2  are  the  measured  autocorrelations  obtained  from  the  observed  data  in 
channels  one  and  two,  respectively,  This  definition  was  first  described  in  Chapter  III. 


K  N 

Ropl  =  if-1  J^41}(x)  l1  ~  5(x)]  +  J>;i(y)]2<Kx)  (4-49) 

k= 1  y=l 

K  N 

RoP2  =  K-1  ^42)(x)[l  -  <J(x)]  +  ^[°p2d(y)]2<5(x)  (4.50) 

1=1  y= 1 

pnew(y)  is  estimated  by  maximizing  the  Q-function  with  respect  to  the  p  param¬ 
eter. 
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dQ 

My0) 


y-y-y-  I  Myp>x)  My0  ~x,x)  p2(y0,x)  ^2(y0  -x,x) 

hhr\  p(y0)  !-p(yo)  1~p( y0) 

-  2o(y0)opi(y0  +  x)  +  2o(y0)op2(y0)  | 

=  0.  (4.51) 

Substituting  in  the  solution  for  onew  from  Eqn.  14.461  and  solving  for  pnew  produces 
a  quadratic  equation  with  the  following  roots: 


pnew(y0) 


B(y0)  +  4/i  +  4^2  ±  ^/[d/i  +  vh2  +  B(y0)]2 

2^(y0) 


45(y0)d>i 


(4.52) 


where 


B(yo)  =  2[S:;i  -  S5]o”(y,).  (4.53) 

The  smallest  (or  positive)  root  is  chosen  when  computing  the  estimate  for  pnew.  The 
sum  of  the  estimated  polarized  objects,  S™{  and  (see  Eqns.  14.441  and  14.45)  are 
unknown  but  easily  computed  by  the  same  maximization  method: 


Qnew 

°opl 

Hnew 

°op2 


(4.54) 

(4.55) 


Finally,  a  solution  is  found  for  onew  using  the  above  expressions  for  pnew, 
and  S™2  and  plugging  them  into  Eqn  14.461  While  this  is  a  complex  equation,  it  is 
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easily  reduced  using  symbolic  mathematical  solver  software  such  as  Mathematica® 
or  MATLAB®.  The  resulting  estimate  for  the  object  is  surprisingly  familiar: 


(4.56) 


Note  the  estimator  is  the  average  of  two  separate  estimates  formed  from  each 
channel.  The  estimator  for  each  channel  is  of  identical  form  of  the  single  channel, 
unpolarized  algorithm  (see  Appendix  [D|.  Essentially,  two  phase  retrieval  estimates 


are  formed  and  averaged  or  fused  together  with  equal  weighting.  Because  autocorre¬ 


lations  are  symmetric  (e.g.  i?/(x)  =  Rj(— x)),  solutions  include  estimates  related  by 
translation  and  180°  rotation.  Therefore,  fusion  of  two  similar  estimates  may  require 
additional  registration  or  alignment  steps. 

4-3  Algorithm  Computation 

In  order  to  initialize  this  iterative  algorithm,  an  initial  guess  for  the  unknown 
object  and  polarization  ratio  is  chosen.  Normally,  the  spatial  bound  of  the  object 
is  known  a  priori  or  estimated  from  the  spatial  extent  of  the  illuminating  beam. 
Also,  the  object  is  positive  and  the  LADAR  system  produces  a  measured  autocorre¬ 
lation.  Therefore,  the  initial  object  guess  is  chosen  with  the  following  conditions:  (1) 
known  support  region  or  spatial  bound  [hi  :  oold(x)  =  0  V  x  9  hi],  (2)  strictly  posi¬ 
tive,  [oold(x)  >  0  V  x  6  fl],  and  (3)  average  value,  A,  computed  from  the  measured 
autocorrelation, 


1/2 


(4.57) 


X=1 


Therefore,  the  initial  object  guess  [oold(x)  :  x  G  fl]  is  selected  as  independent  random 
variables  distributed  uniformly  over  the  interval  [A  —  0.1,  A  +  0.1].  The  initial  guess 
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for  the  polarization  ratio  is  chosen  to  be  uniformly  equal  to  1/2  within  the  support 
region;  [p°ld(x)  =  0.5  V  x  G  fi]. 

4-4  Stopping  Criteria 

With  all  iterative  algorithms,  knowing  when  to  stop  is  an  operationally  com¬ 
pelling  capability.  Often,  algorithms  are  allowed  to  run  for  a  specified  number  of 
iterations  or  as  long  as  operational  time  permits  if  for  each  successive  iteration,  the 
subsequent  estimate  error  is  known  to  be  less  than  the  previous  estimate  error.  In 
some  cases,  noise  amplification  occurs  and  the  estimated  solution  diverges  if  the  al¬ 
gorithm  is  permitted  to  iterate  too  long.  From  simulation  and  experimental  results 
with  exponential  noise,  both  the  two-channel,  polarimetric  phase  retrieval  algorithm 
presented  here  and  Schulz  and  Snyder’s  single  channel  phase  retrieval  algorithm  [38] 
diverge  due  to  noise  amplification  after  too  many  iterations.  Even  with  the  noise  am¬ 
plification,  the  Poisson  model  offers  attractive  properties  as  detailed  in  Refs.  [38, 39] 
and  produces  satisfactory  object  estimates  with  the  optimal  number  of  iterations. 

Throughout  the  literature,  few  algorithms  are  presented  with  stopping  criteria 
based  on  the  statistical  properties  of  the  data.  However,  Phillips  provides  a  simple 
and  compelling  approach  to  this  problem  [33].  Phillips  proposes  a  dampening  routine 
based  on  the  statistics  of  the  data  by  comparing  the  variance  of  the  predicted  object 
to  the  variance  of  the  data  for  each  subsequent  iteration  such  that 

K 

K~v  J^[4(x)  -  R°ld(x)]2  <  (3  •  s2(x),  (4.58) 

k= 1 

where  (3  G  {M  >  1}  is  a  user  chosen  parameter  determining  the  degree  of  dampening 
and  s2  is  the  sample  variance  computed  from  the  observed  data.  If  the  variance  of 
the  predicted  object  is  within  the  variance  of  the  observed  data,  the  pixel  should  be 
damped  by  setting  the  ratio  of  the  measured  object  autocorrelation  to  the  estimated 
object  autocorrelation  equal  to  one, 
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Ro (x)_ 

Rold(*) 


(4.59) 


The  algorithm  is  then  allowed  to  iterate  as  long  as  operationally  feasible.  Phillips 
did  not  propose  a  method  for  selecting  the  optimal  (3  as  this  parameter  is  data  de¬ 
pendent.  When  applying  the  pixel  dampening  criteria  for  this  problem  with  exponen¬ 
tial  noise,  divergence  still  occurred  after  the  optimal  iteration  number  was  exceeded. 
Therefore,  a  global  stopping  criteria  was  selected  to  achieve  the  iteration  number  near 
the  optimal  number  in  terms  of  mean-square  error.  If  the  criteria  found  in  Eqn.  14.581 
is  reached  on  average  throughout  the  entire  data  image,  the  algorithm  is  stopped. 
Selecting  [1.01  <  /3  <  1.3]  produced  satisfactory  results  for  the  simulated  data  set. 

N  K 

(jv  ■!<)-'  £[4(x)  -  iC(x)]2  <  P  ■  («o) 

x=l  k=  1  x 

Summary 

This  chapter  provided  a  detailed  development  for  two  new  phase  retrieval  al¬ 
gorithms  for  the  correlography  problem  presented  in  Chapter  II.  Employing  the  EM 
method  and  using  the  Poisson  noise  model  and  a  polarimetric  data  model  two  iterative 
phase  retrieval  algorithms  were  developed.  Both  MAP  and  ML  methods  were  pre¬ 
sented  as  well  as  computational  considerations.  Finally,  a  statistically  based  stopping 
criteria  was  detailed. 
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V.  Cramer-Rao  Lower  Bound  on  Spatial  Resolution 

This  chapter  provides  a  statistical  analysis  of  the  theoretical  bound  on  resolution 
for  reconstructed  images  computed  in  the  phase  retrieval  problem  presented  in  this 
research.  A  statistical  model  is  presented  along  with  a  simplistic  object  model.  It  is 
often  postulated  how  well  a  particular  approach  performs  as  compared  to  another.  A 
theoretical  bound,  independent  of  the  computational  algorithm  used  is  a  measure  of 
“best”  possible  performance  given  the  underlying  measurements  and  data  statistics. 
The  Cramer-Rao  Lower  Bound  (CRLB)  is  often  used  as  a  such  a  measure.  The  CRLB 
provides  a  lower  bound  on  the  error  covariance  matrix  for  any  unbiased  estimator  [36]. 

Of  interest  with  image  reconstruction  is  how  well  the  image  reconstructed  com¬ 
pares  to  the  original.  An  analysis  of  phase-retrieval  error  is  provided  by  Cederquist 
and  Wackerman  [5],  with  Gaussian  noise  statistics  assumed.  Here  the  variance  of 
the  estimate  of  object  intensity  is  bounded,  not  resolution.  In  addition  to  overall 
reconstruction  error,  resolution  or  differentiation  of  image  detail  is  often  measured  or 
studied.  Resolution  measures  how  well  two  distinct  but  adjacent  objects  or  sources 
may  be  individually  distinguishable  [8,22], 

A  theoretical  bound  on  resolution  for  a  lens-based  imaging  system  was  presented 
by  Shahram  and  Milanfar  [41].  Shahram  and  Milanfar  used  a  two  point  source  model 
for  demonstrating  resolution.  Strong  provided  a  similar  theoretical  bound  for  imaging 
resolution  using  a  two  point  source  model  and  polarimetric  data  [43] .  In  this  research, 
a  similar  approach  is  applied  to  the  non-imaging,  correlography  case.  A  theoretical 
resolution  bound  is  developed  for  three  cases:  (1)  a  single  channel,  unpolarized  system, 
(2)  a  two-channel,  polarimetric  system,  and  (3)  a  dual-channel  polarimetric  system. 
The  second  and  third  cases  will  be  compared  to  the  first  case  in  order  describe  the 
improvement  in  performance  with  the  addition  of  polarization  diversity  in  the  remote 
sensor  system. 
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5.1  Object  Model 

The  object,  o,  consists  of  two  point  sources  separated  by  an  unknown  distance, 
A.  The  viewing  of  this  simple  object’s  autocorrelation  function  is  blurred  by  a  known 
point  spread  function  (PSF)  defined  by  the  observing  aperture.  The  two-dimensional 
object  geometry  considered  here  is  manipulated  in  one  dimension  to  simplify  analysis 
and  computation  as  well  as  maintain  a  well-conditioned  problem.  The  number  of  nui¬ 
sance  parameters  are  minimized.  Adding  additional,  unknown  nuisance  parameters 
(e.g.  spatial  location  and  brightness)  only  increases  the  error  bound  [36].  Thus,  a 
lower  bound  is  computed  when  the  number  of  parameters  is  minimized.  Figure  15.1 
depicts  the  geometry  of  the  object  model. 


v 


Figure  5.1:  Two  Point  Source  Object  Model  Geometry 
The  object  model  is  mathematically  described  by, 

o(u,  v)  =  oi5(u)  +  02<5(m  —  A),  (5.1) 

where  5(u)  is  the  Dirac  delta  function.  For  this  development,  it  is  assumed  the  point 
source  intensities  are  equal,  0\  =  02  with  different  polarization  characteristics.  This 
development  is  easily  extended  to  a  more  complex  case  where  the  point  sources  are 
unknown  and  unequal;  more  unknowns  to  be  estimated  and  a  larger  Fisher  Informa¬ 
tion  matrix.  However,  this  increases  the  number  of  nuisance  parameters  and  moves 
this  towards  an  ill-conditioned  problem;  especially  when  polarization  is  considered. 
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Problem  stabilization,  or  regularization,  must  be  introduced  to  overcome  the  under¬ 
determined  nature  of  the  problem  or  the  Fisher  Information  (FI)  matrix  becomes 
non-invertible.  For  computational  algorithms,  this  is  overcome  in  practice  with  the 
introduction  of  a  priori  information  such  as  object  constraints  (positivity,  spatial 
bound,  etc.)  as  described  in  Chapter  For  the  theoretical  bound  computation 
in  this  development,  the  problem  is  kept  relatively  well-conditioned  by  limiting  the 
number  of  nuisance  parameters  in  the  FI  matrix. 

This  simple,  two-point  object  model  would  not  physically  create  laser  speckle 
at  the  observation  plane  as  a  large  number  of  scatterers  are  required  to  create  the 
speckle  phenomenon.  However,  this  object  model  is  identical  to  the  the  framework 
established  in  the  imaging  case  [41,43].  Also,  this  theoretical  development  employs 
the  statistical  model  described  in  Chapter  III.  This  simple,  two  point  source,  model 
is  repeated  here  for  the  purpose  of  showing  relative  improvement  provided  by  adding 
polarization  diversity.  With  demonstrated  improvement  using  polarization  diversity 
in  the  simple  case,  it  is  postulated  the  improvement  is  similarly  obtained  in  the  case 
of  a  complex  object  with  large  number  of  polarization  diverse  scatterers. 

The  autocorrelation  of  this  object,  blurred  by  a  known  PSF  and  corrupted  with 
exponential  noise,  is  used  to  depict  a  theoretical  resolution  bound  with  and  without 
polarization.  In  the  correlography  case,  the  PSF  is  known  with  no  atmospheric  cor¬ 
ruption  observed  (see  Chapter  11.2.21).  However,  the  PSF  does  limit  the  resolution  of 
the  estimated  object  as  only  a  spatially  limited  and  discrete  sample  of  the  speckle 
intensity  is  collected.  The  detection  array  is  modeled  with  a  square,  uniform  aperture 
with  width,  S,  or  equivalently  a  rectangular  (rect)  function, 


( 


0,  |v|  >  f, 


(5.2) 
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where  v  is  the  two-dimensional  coordinate  vector  in  the  detector  plane.  Without  addi¬ 
tional  corruption  such  as  atmospheric  distortion,  the  PSF  of  the  detector  array,  g(x), 
is  the  Fourier  transform  of  the  array  aperture  function,  A(v),  magnitude  squared, 


<7(x)  =  |4{A(v)}|2 

=  S2 sine2  ,  (5.3) 

where  T  represents  the  two-dimensional  Fourier  Transform  and  x  and  v  are  two- 
dimensional  coordinate  vectors. 

Let  c4(x)  denote  the  kt\i  frame  of  observed  data  (a  single  autocorrelation  cor¬ 
rupted  by  exponential  noise).  The  expected  value,  P0(x),  of  the  observed  data  will 
be  the  autocorrelation  of  the  object  convolved  with  the  PSF,  g : 


i?u(x)  =  /  /  o(u)o(v  +  u)g(x  —  v)dudv, 


—  OO 
OO 


[oi5(u)  +  o25(u  —  A)]  [oi<5(v  +  u)  +  o2<f(v  +  u  —  A)]  g(x  —  v)dudv 
=  (°i  +  °\ Mx)  +  0i02  [  g(x  -  A)  +  g(x  +  A)  ] .  (5.4) 


—  OO 

J2  ,  J2\ 


It  is  assumed  each  noisy  autocorrelation  frame  is  statistically  independent  and 
the  noise  in  each  frame  is  statistically  independent  at  every  sampled  data  point. 
Therefore,  the  joint  probability  density  function  is  formed  and  detailed  as 


p[di(x), ...,  d*(x)] = n  n  exp 

k  x  °V  / 


-4(x) 

-Ro(x) 


(5.5) 


69 


5.2  Resolution  Criterion 


Resolution  criterion  must  be  defined  in  order  to  draw  any  meaningful  conclu¬ 
sions  or  make  reasonable  comparisons.  “Two-point  resolution,  which  is  defined  as  the 
system’s  ability  to  resolve  two  point  sources  of  equal  intensity,  is  a  widely  used  mea¬ 
sure  of  the  overall  resolving  capabilities  of  an  imaging  system  [8].”  Using  statistical 
analysis,  the  two-point  resolution  model  is  adopted  here. 

Let  cr^  denote  the  lower  bound  on  the  mean-squared  error  for  any  unbiased 
estimator  of  A.  Equivalently,  a  a  is  the  lower  bound  on  standard  deviation,  or  root 
mean  square  deviation,  for  an  unbiased  estimator.  For  this  research,  two  closely- 
spaced  point  sources  are  considered  resolved  (e.g.  distinguishable  from  a  single  point) 
if  the  separation,  A,  is  greater  than  one  standard  deviation,  a  a,  of  the  estimate  error: 


A  >  <j a-  (5.6) 

If  the  two  points  are  separated  by  less  than  one  standard  deviation,  the  uncertainty 
of  the  estimate  is  large.  In  this  case,  the  uncertainty  of  the  estimated  separation,  on 
average,  tends  larger  than  or  approximates  the  actual  separation.  If  the  two  points  are 
separated  by  greater  than  one  standard  deviation  of  estimate  error,  the  uncertainty 
of  the  estimate,  on  average,  tends  to  be  much  less  than  the  actual  separation.  This 
criterion  appears  somewhat  arbitrary:  two  standard  deviations  could  reasonably  be 
chosen.  However,  this  criterion  is  selected  for  comparing  correlography  systems  with 
and  without  polarization.  This  criterion  is  identical  to  the  imaging  case  found  in 
Strong  [43]. 

5.3  Bound  for  Single,  Unpolarized  Channel 

First,  a  CRLB  is  computed  for  a  single- channel,  unpolarized  system.  For  com¬ 
puting  the  lower  bound,  consider  the  case  where  o\  =  o2.  The  data,  c4(x),  is  observed 
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without  a  polarizing  reference  or  analyzer.  The  log-likelihood  function,  L,  is  detailed 


as 


l(a,0i)  =  ££ 

k  x 


111  R,0  (x) 


4(x) 

-Ro(x) 


The  FI  matrix,  J,  is  calculated  by 


(5.7) 


*7 


'  92L  ' 
_dLidLj 


where  E[-]  is  the  Expected  Value  operator.  Evaluating  Eqn.  15.81  yields: 


(5.8) 


Jij  —  —E 


d  d 
f)L,  <)L 


EE 


=  -E 


d 

dLi 


-*E 


In  R0(x) 


dR0(x) 


Rn  (x)  dL,, 


4( x) 

-Rn(x) 


EE 


4(x)  9i?Q(x) 
i?2(x)  dLj 


=  -E 


K  \ 

1  dR0(x)  dRa(x) 

1  d2R0{x)' 

L  X 

_i?2(x)  <94  94 

R0(x)  dLidLj  _ 

■EE 

fc  : 

*E 


24 (x)  dR0(x)  dR0(x)  +  4(x)  d2Ra(x 


4KX)  dLi 
— 1  dR0(x)  dR0(x 


dLj 


i?2(x)  94 


9L, 


+ 


-R2(x)  dLidLj 
1  d2RJx) 


i?0(x)  dLidLj 


EE 


'2 R0(x)  dRa(x)  dRa(x) 

R0(x) 

d2R0{x)~ 

_  R3(x)  dLi  dLj 

^o(X) 

dL^Lj  _ 

a'E 


1  dRa(x)  dR0(x) 
i?2(x)  94  <94 


(5.9) 


Evaluating  the  partial  derivatives  of  the  autocorrelation 
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3^X'  =  |2oi9(x)  +  °i  (»(x  -  A)  +  »(x  +  A))  J 


(5.10) 


dRa(x)  d 
do  i  doi 


4ois(x)  +  2 o,  [s(x  -  A)  +  g(x  +  A)] 


(5.11) 


As  stated  previously,  the  bound  with  point  source  intensity  0\  =  02  is  com¬ 
puted.  In  this  case,  one  must  jointly  estimate  two  parameters:  A  and  o This 
produces  a  2  x  2  Fisher  Information  matrix.  Because  this  is  an  ill-posed  problem, 
there  are  multiple  solutions  possible  and  the  Fisher  Information  matrix  may  be  ill- 
conditioned.  Adding  more  than  four  unknowns  causes  the  Fisher  Information  matrix 
to  be  highly  ill-conditioned  and  approach  singularity.  The  math  is  easily  extended 
to  more  unknowns  (e.g.  o\  7^  02,  two-dimensional  separation);  however,  the  addition 
of  nuisance  parameters  when  including  polarization  produces  a  highly  ill-conditioned 
matrix  where  its  inverse  is  not  meaningful.  Throughout  this  development,  to  include 
comparison  to  cases  with  polarization,  the  Fisher  Information  matrix  is  kept  as  small 
as  possible  to  create  a  reasonably  conditioned  matrix.  Computing  the  elements  of 
the  Fisher  Information  matrix  yields  the  following: 
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Jaa  =  kJ2 


-R^(x)  \  <9  A 


dRJx) 


=  K  V  — 1 

x 

Joio! =  k  'y  ^ 


i?o(x)  [dx 
1 


■7-5 (x  +  A)  -  -^-g(x  -  A) 

GtX 


<9i?0(x) 


=  *E 

X 

ja01  =  kJ2 


i?2(x)  V  9oi 
1 


n  2 


4oi0(x)  +  2oi  [#(x  -  A)  +  g(x  +  A)] 


^(x 
1  dR0(x)  dRa(x 
R%  (x)  <9o!  <9  A 


A'Efi2(x){ 


.  ^(x 

2 of  g(x  -  A)  +  s(x  +  A) 


4o?^(x)  [g'(x  +  A)  -  c/'(x  -  A)] 

g'(x  +  A)  -  g'(x  -  A) 


Jq\A  Jaoi 


(5.12) 


(5.13) 

(5.14) 


The  Fisher  Information  matrix  is  constructed  as 


J  = 


JaA  J A01 
Jo\A  Jo\o\ 


(5.15) 


The  lower  bound  on  the  error  Covariance  matrix,  C,  is  found  by  taking  the 
inverse  of  the  Fisher  Information  matrix: 


C  >  J1, 


(5.16) 


and  the  CRLB  for  the  separation  parameter,  A,  is  [36] 
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<4> 


r— l 


(1,1) 


The  two  point,  sources  are  considered  resolved  if 


(5.17) 


A  > 


(5.18) 


This  simple  case  is  considered  the  baseline  and  will  be  used  comparatively  with  po¬ 
larimetric  cases. 


5-4  Bound  for  Two-Channel,  Polarimetric  Estimator 

Next,  the  resolution  bound  for  a  two-channel  polarimetric  estimator  is  calcu¬ 
lated.  A  two-channel  system  observes  both  polarized  and  unpolarized  data  in  two 
independent  channels.  If  the  data  observed  in  the  polarized  and  unpolarized  channels 
are  statistically  independent,  the  joint  probability  density  function  can  be  expressed 
as  a  product  of  the  marginal  density  functions.  The  joint  PDF  is 


p[4(x),d;(x)] = nnn 

k  l  x 


1 

P0i(x)Po2(x)  CXP 


-4(x)  di  (x) 
i?ol(x)  Rq2  (x) 


(5.19) 


where  R0\  is  the  expected  value  of  the  data  observed  in  the  unpolarized  channel  and 
R0 2  is  the  expected  value  of  the  data  observed  in  the  polarized  channel.  The  polarized 
channel  observes  the  autocorrelation  of  the  object  as  viewed  through  the  polarizer. 
The  polarized  object,  op,  is  described  by 
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Op(u)  =  p(u)o(u) 


=  oipi5(u)  +  o2p2<S(u  -  A).  (5.20) 

The  expected  value  of  the  observed  data  in  the  unpolarized  channel  remains  as  de¬ 
tailed  in  Eqn.  15.41  The  expected  value  of  the  observed  data  in  the  polarized  channel 
is 


R02 (x)  =  //  Op(u)op(v  +  u)p(x  -  v)dudv 


—  OO 

oo 


Oipi<5(u)  +  o2p2<5(u  -  A) 
0iPi<5(u  +  v)  +  o2p2<5(u  +  v  -  A) 


g(x  —  v)dudv 


oipfS(u)S(u  +  v)  +  oipio2p2<5(u)5(u  +  v  -  A) 


+  oipio2p25(u  —  A)<5(u  +  v)  +  o2p25(u  —  A)5(u  +  v  —  A) 
°iPi^(v)  +  oipio2p25(v  -  A)  +  0ipi02p2^(v  +  A) 

+  o22pl5(v)  g(x  —  v)dv 

=  [o\p\  +  o22p22]g{x.)  +  0\p\02p2  [p(x  -  A)  +  g(x  +  A)] . 


g(x  —  v)dudv 


(5.21) 


If  Oi  =  o2  then  the  expected  values  in  the  unpolarized  and  polarized  channels  simplify 
to 
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^oi(x)  =  2ofo(x)  +  o\  [g(x  -  A)  +  g{x  +  A)] 

Ro 2(x)  =  o\  [p\  +  pi]  g(x)  +  olp!p2  [s(x  -  A)  +  g(x  +  A)] . 


(5.22) 

(5.23) 


The  log-likelihood  function,  L  is  detailed  as 


L(A,o1,p1,p2)  =  ZEE  —  In  R0\  (x)  —  In  R, 

hi  T.  ^ 


'o2  X  - 


4(x) 


d;(x)  \ 

i?o2(x)  J  ' 


(5.24) 


Computing  the  Fisher  Information  matrix  yields 


Jij  —  —E 


d2L(A1oi,pi,p2) 


dLidLj 


=  K2 


X 


1  <TRo:l(x)  dR0i{x) 


Rl 


oil 


DU  ()L, 


1  dRo2(x)  dRo2(yi 
R2o2(x)  dLi  dLj 


(5.25) 


Evaluating  the  partial  derivatives  of  the  autocorrelation  in  the  polarized  channel 


yields 


dR0 a(x) 

dA 


dR02{y) 

do\ 

dR02{y) 

dpi 

dR02  iy) 

dp2 


[pi  +  pI\  #(x)  +  °iPiP2  [#(x  -  A)  +  g(x  +  A)] 


o\piP2  </(x  +  A)  -  #'(x  -  A) 


2oi  [Pi  +  P2]  0(x)  +  2o1p1p2  [s(x  -  A)  +  c/(x  +  A)] , 
2o?pi^(x)  +  o2p2  [#(x  -  A)  +  g(x  +  A)] , 

2o\p2g(yP)  +  o\pi  [#(x  -  A)  +  g(x  +  A)] . 


(5.26) 

(5.27) 

(5.28) 

(5.29) 
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For  comparison  to  the  single-channel,  unpolarized  case  the  bound  with  point 
source  intensity  Oi  =  o2  is  computed.  The  two  point  sources  may  not  have  the 
same  degree  of  polarization  and  will  produce  different  intensities  as  viewed  through  a 
polarizer.  The  unknown  parameters  consist  of  the  separation  parameter,  A,  the  point 
source  intensity,  o1;  and  the  polarization  ratio,  p.  Again,  the  number  of  nuisance 
parameters  is  minimized  to  avoid  an  ill-conditioned  matrix.  The  elements  of  the 
Fisher  Information  matrix  are  calculated  next. 
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Jaa  =  K2J2 
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01P1P2 
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(5.30) 


&R0i(x)>\2+  1  fdRo2(x)'2 
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<9c>i 


4oip(x)  +  2oi 


^02  (x) 


<9c>i 
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+a-e 
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2oi  [p2  +  P2]  fi'(x)  +  MP1P2 


^(x  -  A)  +  g(x  +  A) 


(5.31) 
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^02  (x) 


2o2pi5f(x)  +  o2p2 
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(5.40) 

(5.41) 
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(5.42) 

(5.43) 

(5.44) 

(5.45) 

(5.46) 


As  stated  previously,  the  error  bound  or  smallest  possible  variances  of  the  es¬ 
timated  parameters  are  calculated  by  taking  the  inverse  of  the  Fisher  Information 
matrix.  The  bound  for  each  unknown  parameter  is  found  along  the  diagonal  of  the 
inverted  matrix.  The  lower  bound  (variance)  on  the  separation  parameter  A,  is  de¬ 
tailed  by 


(5.47) 


5.5  Bound  for  Dual- Channel,  Polarimetric  Estimator 

Next,  the  resolution  bound  for  a  dual-channel  polarimetric  estimator  is  calcu¬ 
lated.  A  dual-channel  system  observes  polarized  and  data  in  two  independent  chan¬ 
nels.  The  polarizer  in  channel  one  is  orthogonal  to  the  polarizer  in  channel  two.  If  the 
data  observed  in  the  two  channels  are  statistically  independent,  the  joint  probability 
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density  function  is  identical  to  Eqn.  15.19  with  the  expected  values  of  the  two  channels 
formulated  differently.  The  joint  PDF  is 


p[4(x),^(x)] = nnn 

k  l  x 


1 

P0i(x)Po2(x)  °XP 


-4(x) 
Pol  (x) 


~di(x) 

Rq2  (x) 


(5.48) 


where  Pol  is  the  expected  value  of  the  polarized  data  observed  in  the  channel  one 
and  R02  is  the  expected  value  of  the  polarized  data  observed  in  the  second  channel. 
With  the  dual-channel  system,  the  two  channels  of  observed  data  are  related  by  the 
following  equation: 


p2  =  l  -Pi,  (5-49) 

where  p\  represents  the  polarization  ratio  of  channel  one  and  p2  represents  the  po¬ 
larization  ratio  of  channel  two.  In  each  channel,  a  noisy  object  autocorrelation  is 
observed  as  viewed  through  the  polarizer.  The  polarized  object  as  viewed  in  channels 
one  and  two  respectively  are  described  by 


Opi(u)  =  pi(u)o(u),  (5.50) 

op2(u)  =  p2(u)o(u).  (5.51) 

Substituting  Eqn.  15. Hand  Eqn.  15.491  into  the  above  equations  for  the  polarized  object 
yields 
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Opi (u)  =  0iPuS(u)  +  o2pnS(u  -  A), 
Op2(u)  =  OiP2l5(u)  +  o2p22<5(u  -  A), 


(5.52) 


=  01  [1  -pn]<5(u)  +  02 [1  —  Pi2] A(u  -  A), 


(5.53) 


where  the  subscripts  on  the  polarization  ratio  (n)  and  (12)  represent  channel  one  and 
object  point  source  one  and  two  respectively.  In  order  to  harmonize  with  previous 
developments,  it  is  assumed  the  point  source  intensities  are  identical,  01  =  o2;  how¬ 
ever,  each  point  source  may  have  a  different  polarization  response,  pu  7^  p\2.  With 
the  polarized  object  described,  the  expected  value  of  the  observed  data  in  the  two 
channels  are  related  to  the  polarized  object  autocorrelation: 


Roi  (x)  =  ol  [p2n  +  p\2\  g(x)  +  olpupu  [p(x  -  A)  +  g(x  +  A)] , 
i?o2(x)  =  o\[(l  -Pu)2  +  (1  -p12)2]c/(x) 


(5.54) 


+  o\{  1  -  pn)(l  -  p12)  [^(x  -  A)  +  g(x  +  A)] . 


(5.55) 


The  log-likelihood  function  and  the  Fisher  Information  matrix  is  constructed 
identical  to  Eqns.  15.241  and  5.251  Evaluating  the  partial  derivatives  of  the  polarized 
autocorrelations  is  similar  to  the  polarized  channel  found  in  two-channel  case.  The 
elements  of  the  Fisher  Information  matrix  are  calculated  next. 
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(5.57) 
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The  lower  bound  on  the  separation  parameter  is  computed  identical  to  the  previous 
cases,  inverting  the  above  Fisher  Information  matrix. 

5.6  Bound  Results  and  Comparison 

For  this  problem,  there  are  challenges  associated  with  computing  the  Fisher 
Information  matrix  and  its  inverse.  First,  from  the  exponential  distribution,  division 
by  the  mean  squared  term  produces  numeric  challenges.  The  sine-squared  PSF  in¬ 
troduces  extreme  nulls  as  it  tapers  to  zero  causing  division  by  very  small  numbers. 
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As  the  separation  parameter  is  varied,  the  nulls  of  the  mean  squared  term  vary.  The 
numerator  also  tapers  to  zero  in  a  similar  manner;  however,  without  modification, 
the  resulting  division  produces  erratic  oscillations  in  the  computed  variance  as  A  is 
varied.  This  challenge  was  eliminated  by  a  computational  mask.  FI  matrix  values 
where  division  by  very  a  small  number  occurs  (below  a  threshold)  are  set  to  zero. 
This  computational  mask  eliminated  the  oscillations  in  the  data  allowing  for  the  vari¬ 
ance  to  decrease  without  discontinuity  as  A  increases.  The  computational  mask  was 
identical  for  all  scenarios  allowing  comparisons. 

Second,  there  are  a  myriad  of  possible  parameter  conditions  producing  varying 
results.  The  computational  scenarios  are  minimized  by  keeping  select  values  constant 
throughout  the  bound  computations.  For  various  bound  computations,  the  parameter 
values  shown  in  Table  [STB' were  chosen  to  aid  in  computation,  comparison  and  analysis. 
For  simplicity,  the  polarization  ratio,  p±,  of  the  first  point  source  is  chosen  to  be  aligned 
with  the  polarizer  (e.g.  p\  =  1).  For  the  dual-channel  case,  pn  is  aligned  with  the 
polarizer  in  channel  one  and  orthogonal  to  channel  two.  The  polarization  ratio  of  the 
second  point  source,  p2  (or  pi2)  is  varied  to  allow  for  several  diversity  scenarios. 

Table  5.1:  Parameter  Values  for  Bound  Computation 


Parameter 

Value 

Ol 

1 

02 

1 

Pi 

1 

P2 

0.25,  0.5,  0.75 

K  Number  of  Frames 

100 

Matrix  Size  (pixels) 

512  x  512 

Aperture  Size  (pixels) 

128  x  128 

The  units  for  the  various  parameter  values  are  briefly  considered.  The  object 
parameter,  o,  is  considered  a  brightness  or  intensity  normally  described  in  photon 
counts  when  measured  with  optical  detectors.  However,  the  photo  count  unit  is 
dropped  due  to  the  transformation  from  laser  speckle  images  to  an  autocorrelation 
image.  The  transformation  changes  the  positive  integer  data  set,  N  to  positive  rational 
numbers,  Q+.  The  estimated  object  strength  is  related  to  the  original  object  via  a 


87 


scale  factor  that  is  a  function  of  LADAR  range  equation  (laser  power,  distance,  etc.), 
detector  integration  time,  aperture  size,  etc.  Because  of  the  transformation,  photo 
count  is  dropped  and  (scaled)  intensity  is  considered.  In  any  experimental  detection 
scenario,  sufficient  photo  counts  must  be  detected  to  produce  fully  formed  speckle  with 
the  appropriate  distribution  related  to  the  illuminated  object  surface.  This  is  assumed 
to  be  properly  accounted  in  system  design.  Extremely  low  light  levels  or  partially 
formed  speckle  conditions  are  not  considered.  The  polarization  ratio  is,  of  course, 
unitless  since  it  is  formed  by  a  ratio  of  polarized  and  unpolarized  object  intensities. 
The  aperture  size  is  detailed  in  number  of  pixels.  Sample  size  may  be  computed  from 
matrix  size  and  detection  geometry.  For  resolution  bound  computations,  matrix  and 
aperture  size  is  selected  without  geometry  consideration.  It  is  assumed  the  sample 
size  is  sufficient  to  meet  critical  sampling  requirements  for  applicable  geometries. 

The  last  computation  consideration  is  a  separation  (or  matrix  shift)  value  of 
less  than  one  pixel.  This  was  accomplished  using  a  MATLAB®  subroutine  that 
implemented  a  sub-pixel  shift  performed  with  the  two-dimensional  Digital  Fourier 
Transform.  The  routine  implemented  a  digital  representation  of  the  Fourier  Shift 
Theorem  [18].  A  phase  shift  of  less  than  one,  (0  <  a  <  1),  was  introduced  in  the 
Fourier  domain  producing  a  sub-pixel  shift  in  spatial  domain.  This  is  essentially 
restated  as 


g(x  -  a:y  -  6)  =  T  1  |G(/*,  fy)  exp[-j2vr(/xa  +  fyb)\  j .  (5.73) 

5.6.1  Single- Channel,  Unpolarized  System.  Figure  15.21  depicts  the  com¬ 
puted  lower  bound  versus  separation  for  the  single  channel  case  using  the  exponential 
statistics  model.  The  aperture  size  is  maintained  constant  and  only  the  separation 
parameter  is  varied.  This  is  the  baseline  for  comparison  with  polarization  diversity 


cases. 


Separation,  A  (pixels) 


Figure  5.2:  Resolution  Bound  vs.  Separation  for  Single  Channel  Case 


Correlography  and  phase  retrieval  performance  does  depend  upon  observing 
aperture  size.  A  larger  aperture  (and  all  other  factors  constant)  provides  smaller 
sample  size  directly  impacting  resolution  performance.  Figure  15.31  demonstrates  the 
effect  of  changing  the  aperture  size  (in  pixels).  As  depicted,  the  bound  computation 
performs  as  expected. 


Figure  5.3:  Resolution  Bound  vs.  Separation  for  Various  Aperture  Sizes 

5.6.2  Two-Channel,  Polarimetric  System.  Next,  the  resolution  bound  for 
the  two-channel  system  is  computed  and  compared  with  the  single  channel  system. 
First,  it  will  be  assumed  the  unpolarized  channel  in  the  two-channel  system  has  the 
identical  light  level  compared  to  what  is  collected  in  the  single-channel  system.  Fig¬ 
ure  15.41  depicts  the  computed  bound  for  the  two  cases  as  a  function  of  the  separation 
parameter,  A.  Note,  the  improvement  demonstrated  in  the  two-channel  case  is  ex¬ 
actly  10  times  the  bound  of  the  single  channel  case.  Because  the  computation  was 
performed  with  K  =  100  frames,  the  improvement  is  exactly  \/~K.  This  improvement 
is  directly  related  to  the  independence  of  the  the  two  channels  and  the  multiplicative 
1/A"  factor  in  the  variance  computation.  No  significant  change  in  performance  was 
observed  by  changing  the  polarization  difference  between  the  two  point  sources  in  the 
two  channel  system.  Therefore,  the  conclusion  drawn  from  this  result  is  the  secondary 
channel  provides  significant  improvement  if  statistical  independence  is  observed.  It 
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can  be  surmised  that  partial  dependence  or  correlation  between  the  two  channels 
impedes  the  improvement  and  may  be  a  subject  of  future  research. 


Figure  5.4:  Resolution  Bound  vs.  Separation  for  One  and  Two  Channel  Cases 

Only  minute  differences  (~  ICO7)  are  observed  by  changing  the  polarization 
ratio  of  the  second  point.  The  minute  changes  do  demonstrate  slight  improvement 
when  more  light  is  observed  (less  suppression  due  to  polarimeter  effects).  Due  to 
only  very  small  computational  differences,  it  is  concluded  polarization  differences  in 
the  object  scene  are  not  important  for  this  simple  two-point  object  model  and  the 
two-channel  system. 

As  detailed  in  Section  [3 .1.11  S  +  P  results  in  a  data  set  not  statistically  indepen¬ 
dent  from  either  S  or  P  channel  data.  Therefore,  the  same  bound  computation  was 
repeated  for  the  two-channel  case  with  one-half  the  light  level  representative  of  collec¬ 
tion  with  a  non-polarizing  beam  splitter  and  polarization  analyzer.  The  polarization 
difference  in  the  two  points  was  maintained  the  same.  The  resulting  bound  curve 
was  found  to  be  identical  to  Figure  E3j  as  expected  with  an  exponential  noise  model. 
The  SNR  for  an  exponential  noise  model  is  exactly  one  (mean  equal  to  standard  de- 
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viation);  therefore,  resulting  bound  computations  are  expected  to  be  independent  of 
signal  strength. 

5.6.3  Dual-Channel,  Polarimetric  System.  Figure  [S3]  depicts  the  resolution 
bound  for  the  dual-channel  system  as  a  function  of  the  separation  parameter.  Multiple 
curves  are  shown  with  three  different  cases  of  polarization  difference  between  the  point 
sources.  Of  important  note,  the  closer  the  two  point  sources  are  in  polarization  ratio, 
the  more  light  is  collected  and  the  bound  is  improved  (lower).  It  is  clear  from  this 
result,  the  polarization  diversity  scheme  depends  mostly  upon  statistical  independence 
in  the  data  channels  and  less  so  on  polarization  difference  between  the  point  sources. 
This  result  is  specific  to  the  two  point  source  model;  however,  generalized  conclusions 
can  be  surmised.  For  complex  objects,  statistical  independence  in  the  data  channels 
is  of  primary  importance  in  designing  a  polarimetric  correlography  system.  Also, 
detection  schemes  should  be  designed  to  maximize  light  levels.  For  example,  a  PBS 
is  superior  to  a  standard,  non-polarizing  beam  splitter  with  polarization  analyzer  due 
to  transmitted  light  levels. 


Separation,  A  (pixels) 

Figure  5.5:  Resolution  Bound  vs.  Separation  for  Dual  Channel  with  Various  Po¬ 
larization  States 
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Figure  15.6  compares  the  performance  of  the  two-channel  system  to  the  dual¬ 
channel  system  for  the  same  polarization  difference.  The  two-channel  system  demon¬ 
strates  slightly  better  performance;  however,  this  can  be  attributed  to  less  amplitude 
suppression  due  to  the  polarimeter  effects  in  the  observed  data. 


Separation,  A  (pixels) 


Figure  5.6:  Resolution  Bound  vs.  Separation  for  Two  Channel  and  Dual  Channel 
Cases  (pn  =  l,pi2  =  0.5) 


The  bound  computations  demonstrate  both  polarimetric  systems  outperform 
the  unpolarized  system  by  a  factor  approximately  equal  to  \/~K.  This  is  primarily 
due  to  statistical  independence  of  the  two  observed  data  channels.  Also,  the  bound 
computations  for  the  simple,  two-point  model  demonstrate  scene  diversity  (polariza¬ 
tion  difference)  is  not  an  important  factor  for  reducing  the  variance  of  an  unbiased 
estimator.  The  phase  retrieval  problem  is  inherently  ill-posed  and  estimators  designed 
to  solve  it  may  be  biased.  Quantifying  bias  remains  an  important  area  of  study  in 
the  field  of  phase  retrieval. 
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VI.  Results  and  Analysis 

This  chapter  provides  simulation  and  experimental  results  generated  using  the  algo¬ 
rithms  detailed  in  Chapter  IHVl  An  analysis  of  the  results  is  provided  with  comparison 
to  a  previously  published  algorithm  [38]. 

6.1  Simulation  Results 

This  section  details  the  results  of  testing  the  Polarimetric  EM  Phase  Retrieval 
algorithm  using  simulated  data.  The  simulated  data  consists  of  100  frames  per  channel 
of  the  object  autocorrelations  corrupted  with  statistically  independent  exponential 
noise  (see  Appendix  A).  With  the  data  formed  in  this  manner,  the  signal-to- noise 
ratio  (SNR)  for  each  individual  frame  is  approximately  one.  The  simulated  object 
was  stored  in  a  128  x  128  discrete  array  with  the  object  embedded  within  the  support 
region  defined  as  the  central  64  x  64  portion  of  the  array. 

The  simulated  object  consisted  of  three  horizontal  bars  of  different  strength  with 
each  bar  assigned  a  distinct  polarization  angle:  bar  one  is  zero  degrees  (6\  =  0°),  bar 
two  is  45  degrees  (#2  =  45°),  and  bar  three  is  60  degrees  (#3  =  60°).  The  distinct  angle, 
9,  represents  fully  polarized  light  oriented  along  the  described  angle  measured  from 
the  vertical  axis.  Figure  16,1  depicts  the  three  bars  with  the  described  polarization 
angles. 
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Bar  1 


Bar  2 


Bar  3 


Figure  6.1:  Simulated  Object  Polarization:  3  Bars  with  Various  Polarization  Angles 


The  object  autocorrelation  was  formed  and  then  subjected  to  a  random  expo¬ 
nential  number  generator  for  100  statistically  independent  realizations  of  the  noisy 
autocorrelation.  The  data  for  the  polarization  channel  was  formed  identically  to  the 
unpolarized  channel  after  computing  the  polarized  version  of  the  object,  op  =  p  x  o. 
The  simulated  polarization  ratio  was  computed  using  the  following  equation: 


p  =  cos2  ((f)  —  6) ,  (6.1) 

where  (f>  is  the  orientation  angle  of  the  polarizer  and  6  is  the  orientation  angle  of  the 
object  element.  Example  data,  with  the  polarizer  transmission  angle  set  to  90°,  is 
pictured  in  Fig.  16.21 
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(C)  (d) 


Figure  6.2:  Example  Simulation  Data:  True  Object  (a),(c)  and  Average  Autocor¬ 
relation  (b),(d)  (K=100)  for  unpolarized  channel  and  polarized  channel  [90°] 

The  Polarimetric  EM  algorithm  and  Schulz  and  Snyder’s  single-channel,  ML 
phase  retrieval  algorithm  [38]  were  both  run  under  identical  conditions  for  comparison. 
The  only  difference  being  the  addition  of  the  polarized  channel  data  and  the  initial 
guess  for  the  p  matrix  needed  for  the  Polarimetric  EM  algorithm.  The  two  algorithms 
were  run  for  50  independent  trials  with  randomly  drawn  initial  guesses  as  described 
in  Section  14731  Total  image  mean-squared  error  (MSE)  [36]  of  the  estimate,  o,  is  the 
metric  of  choice  for  comparing  simulated  performance: 


MSEio)  =£[(o-o)2] 

-  N-1X}[0(x)  -o(x)]2,  (6.2) 

X 


where  N  is  the  total  number  of  pixels. 
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The  iterative  algorithm  estimate  is  dependent  upon  the  initial  guess.  Comparing 
with  identical  iteration  number,  a  slightly  different  estimate  is  provided  for  every 
different  initial  guess.  Because  the  autocorrelation  is  symmetric  and  bounded  (image 
goes  to  zero  at  boundary),  the  allowed  solution  set  includes  estimates  related  by 
translation  and  180°  rotation.  Therefore,  registration  of  the  recovered  object  with 
the  original  object  was  performed  before  computing  error  measures.  This  effort  used 
a  vector-based,  energy-normalized,  non-circular  cross  correlation  technique  for  image 
registration  [2] .  The  error  results  from  the  50  trials  were  averaged  to  form  a  final  MSE 
result.  Each  algorithm  was  allowed  to  run  for  a  specified  number  of  iterations  before 
stopping.  The  polarimetric  EM  algorithm  was  also  run  with  the  global  stopping 
criteria  described  in  Section  14,41 

Figure  16.3  depicts  the  results  of  the  MSE  versus  iteration  number  for  both 
algorithms  with  the  polarizer  transmission  axis  set  to  90°.  The  EM  polarimetric 
algorithm  converges  faster  and  to  a  smaller  MSE  as  compared  to  the  ML  single 
channel  algorithm.  Figure  6.3  also  depicts  the  need  for  the  stopping  criteria.  By 
adding  polarization  sensitivity  to  the  sensor  array,  improvement  in  excess  of  12  percent 
is  observed.  Similar  results  were  observed  with  the  polarizer  transmission  axis  rotated 
to  other  angles. 
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Figure  6.3:  Simulated  Error  Results  for  Two  Algorithms  per  Iteration  Number 

(K=100  frames) 


Figure  16.4  shows  performance  versus  number  of  laser  speckle  frames,  K ,  related 
to  operational  collection  time.  For  the  comparison  in  Fig.  I6.4L  the  ML  single-channel 
algorithm  was  stopped  after  exactly  150  iterations  and  the  polarimetric  EM  algorithm 
stopped  prior  to  150  iterations  near  the  optimum  with  (3  =  1.025.  The  stopping 
mechanism  can  successfully  stop  the  algorithm  near  the  optimum  iteration  number. 
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Figure  6.4:  Simulated  Error  Results  for  Two  Algorithms  per  Frame  Number,  K 

The  dual-channel  algorithm  is  very  similar  but  both  observation  channels  are 
polarized.  Both  channels  provide  polarized  data  but  the  transmission  axis  of  the  two 
channels  are  orthogonal  (90°  rotation).  The  dual-channel  algorithm  was  also  tested 
with  simulated  data  and  compared  to  the  single-channel,  ML  algorithm  and  the  two- 
channel  EM  algorithm.  Figure  [6751  depicts  the  simulated  results  as  a  function  of  iter¬ 
ation  number.  For  the  simple,  3  bar  object,  the  dual-channel  algorithm  also  provided 
improvement  over  the  single-channel  non-polarized  ML  algorithm.  The  dual-channel 
algorithm  estimated  the  two-dimensional  object  with  less  MSE  and  fewer  iterations. 
The  stopping  criteria  was  also  successfully  implemented  with  the  dual  channel  al¬ 
gorithm.  With  proper  selection  of  the  dampening  parameter,  f3,  the  algorithm  can 
be  stopped  near  the  optimum  iteration  number  or  prior  to  the  fixed  150  iterations. 
Figure  [6761  depicts  the  simulated  error  results  as  a  function  of  frame  number,  K.  Com¬ 
parison  to  the  single-channel,  ML  algorithm  was  accomplished  with  a  fixed  iteration 
number  of  150.  The  stopping  mechanism  is  very  sensitive  to  the  selection  of  (3  and 
further  analysis  is  needed  to  determine  a  method  for  selecting  the  optimum  (3. 

The  dual-channel  algorithm  averages  the  estimator  results  for  each  channel  (see 
Eqn.  14.56).  This  fusion  of  the  two  separate  data  sets  can  be  problematic  if  the  data 
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sets  are  not  aligned  or  registered.  A  simple  image  registration  routine  was  employed 
prior  to  adding  the  estimates  in  the  two  channels.  Further  improvement  may  be 
obtained  with  better  data  fusion  techniques. 


Figure  6.5:  Simulated  Error  Results  for  Three  Algorithms  per  Iteration  Number 

(K=100  frames) 
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Figure  6.6:  Simulated  Error  Results  for  Three  Algorithms  per  Frame  Number,  K 


6.2  Experimental  Results 

In  order  to  validate  the  theoretical  work  and  computer  simulations,  a  simple 
laser  speckle  experiment  was  performed.  The  polarimetric  phase  retrieval  algorithm 
successfully  recovered  a  three  bar  object  from  a  series  of  noisy  autocorrelations  formed 
from  collected  laser  speckle  images. 

The  laboratory  experiment  was  conducted  with  available  laboratory  hardware. 
The  laboratory  experiment  is  not  completely  representative  of  the  (proposed)  large- 
scale  system  designed  to  recover  remote  satellite  images;  however,  it  does  serve  to 
reinforce  the  theory  and  development  presented  in  this  research.  The  choice  of  labo¬ 
ratory  hardware  was  purely  out  of  convenience  and  availability.  Much  improvement 
in  experimental  performance  is  available  via  hardware  and  experimental  design.  Even 
so,  the  chosen  laboratory  hardware  does  perform  reasonably  similar  in  function  to  the 
proposed  sensor  system. 

A  traditional  charge-coupled  device  (CCD)  camera  and  polarizing  film  was  used 
to  observe  two  channels  of  data  from  a  coherently  illuminated  object.  For  this  lab- 
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oratory  experiment,  the  hardware  configuration  is  simplified  by  using  a  continuous 
wave  (CW)  source  vice  a  pulsed  laser  source.  Sufficient  light  level  is  achieved  by 
appropriate  integration  time  at  the  camera. 


The  laboratory  setup  is  depicted  in  Fig.  16.71  The  test  sensor  consists  of  a  Pho¬ 
tometries  Cascade  512B  camera  without  a  lens,  with  removable  polarization  analyzer 
placed  in  front  of  the  camera  aperture.  The  camera  is  an  electron-multiplying  charge- 
coupled  device  sensor.  The  camera  array  is  512  x  512  pixels  with  a  16/irn  pitch.  A 
spatially  coherent  source  at  630  nm  was  used  to  back  illuminate  a  target  set.  A  laser 
line  filter  at  630  nm  was  inserted  at  the  camera  aperture  to  minimize  background 
light. 
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Figure  6.7:  Diagram  of  Laboratory  Setup 


The  target  set  consists  of  a  glass  plate  completely  opaque  where  no  object 
exists  and  transparent  where  the  object  exists.  The  experimental  object  consisted 
of  three  identically  sized  bars.  To  emulate  polarization  effects,  the  side  opposite 
illumination  of  two  of  the  bars  were  covered  with  polarizing  film  aligned  to  the  same 
axis.  To  emulate  random  surface  roughness,  a  randomizing  phase  screen  was  placed 
between  the  source  and  the  bar  target  consisting  of  highly  fibrous,  white  paper.  The 
paper  phase  screen  was  moved  for  each  collection  frame  to  emulate  random  phase 
perturbations  and  produce  statistically  independent  laser  speckle  images. 

The  propagation  path  included  a  0.5m  lens  at  the  target  plane  to  emulate  far- 
field  (or  Fraunhofer)  propagation  to  the  camera  array  plane.  Laser  power  was  not 
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a  consideration  as  the  source  was  placed  near  the  target  plane.  Camera  integration 
time  was  selected  after  several  preliminary  tests  and  then  held  constant  throughout 
the  experiment. 

For  extremely  remote  sensing  such  as  imaging  of  space-borne  objects,  the  over¬ 
all  path  length  is  very  large  compared  to  distance  where  atmosphere  turbulence  is 
encountered  enabling  us  to  consider  the  atmospheric  turbulence  as  a  single,  uniformly 
distributed  phase  screen.  In  this  research,  a  layered  atmosphere  and  scintillation  are 
not  considered.  Experimentally,  a  small  path  length  is  used  to  ensure  only  uniform 
atmosphere  is  encountered. 

In  an  operational  system  with  LADAR  backscatter,  the  target  geometry,  sur¬ 
face  roughness  and  propagation  distances  typically  produce  the  desired  laser  speckle 
effects.  For  the  laboratory  experiment,  the  paper  phase  screen  enabled  experimen¬ 
tation  using  back  illumination  and  simplified  the  overall  experiment.  See  Ref.  [44] 
for  experimentation  with  a  backscatter  setup.  The  paper  phase  screen  produced  low 
light-levels  where  read  noise  dominates  the  detection  process,  though  sufficient  for  the 
experiment.  However,  the  paper  did  exhibit  spatially  dependent  surface  roughness 
for  the  spatial  sampling  size  produced  by  the  camera.  The  correlography  technique 
assumes  spatially  independent  surface  roughness  (see  Sec.  II. 2. 111.  To  overcome  the 
effects  of  spatially  dependent  surface  roughness  and  low-light,  levels,  each  512  x  512 
laser  speckle  image  was  segmented  into  16,  128  x  128  laser  speckle  images.  Each 
128  x  128  laser  speckle  image  contains  the  statistical  nature  of  the  target  and  can  be 
processed  independently.  Target  spatial  resolution  is  lost  but  SNR  is  improved  by  a 
factor  of  four  and  spatially  independent  surface  roughness  is  gained  due  to  coarser 
sampling  in  the  target  plane.  In  order  to  minimize  reflections  from  lab  equipment 
and  background  light,  the  propagation  path  was  enclosed  in  a  light  baffle.  An  image 
of  the  back  illuminated  target  taken  with  the  test  camera  (with  lens  and  without 
polarizer)  of  the  bar  target  is  shown  in  Fig.  6.8[  The  image  depicts  the  bright  center 
bar  and  the  two  side  bars  with  reduced  brightness  due  to  the  effects  of  the  polarizing 
film. 
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Figure  6.8:  Image  of  Back  Illuminated  Bar  Target  Set 

Figure  6.9  depicts  the  recovered  image  using  the  polarimetric  phase  retrieval 
algorithm.  The  recovered  image  was  produced  after  28  iterations.  The  three  bar 
target  set  is  clearly  depicted  with  the  side  bars  reduced  in  intensity  compared  to 
the  center  bar.  Figure  16.101  depicts  the  recovered  image  using  the  single-channel 
algorithm  developed  by  Schulz  and  Snyder  [38].  For  visual  comparison,  this  image 
was  28  iterations  and  using  the  same  starting  guess  as  the  two-channel  solution.  For 
the  single-channel,  Schulz-Snyder  algorithm  and  the  same  iteration  number,  three 
bars  are  clearly  discernable;  however,  the  bar  shape  is  more  rounded  and  less  defined 
compared  to  the  two-channel  polarimetric  algorithm. 
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Figure  6.9:  Recovered  Image  From  Experimental  Data;  Two-Channel  Algorithm, 

28  Iterations 


Figure  6.10:  Recovered  Image  From  Experimental  Data;  Single- Channel  Algorithm, 
28  Iterations 
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VII.  Conclusion 


Previous  chapters  presented  a  detailed  background  and  motivation  for  correlography 
and  phase  retrieval  that  may  be  employed  in  a  specific  remote  sensing  scenario.  This 
research  is  primarily  motivated  by  the  need  to  produce  two-dimensional  images  of 
remote  satellites  from  earth-bound  sensors.  For  satellites  located  at  geosynchronous 
orbit,  traditional  imaging  schemes  are  inadequate  to  overcome  distance  and  atmo¬ 
spheric  distortion.  Additionally,  remote  sensing  of  the  earth’s  surface  from  space- 
borne  systems  may  also  be  improved.  The  primary  motivation  for  considering  pupil 
plane  imaging  is  that  very  large  apertures  can  be  synthesized  with  a  large  collection 
of  simple,  inexpensive  detector  elements  without  a  monolithic  lens.  The  correlog¬ 
raphy  technique  using  pupil  plane  imaging  coupled  with  appropriate  phase  retrieval 
algorithms  has  shown  potential  to  provide  cost  effective  systems  for  defeating  dis¬ 
tance  and  atmospheric  hurdles.  This  research  provided  a  new  investigative  initiative 
to  solve  this  difficult  problem.  This  chapter  concludes  the  main  document  by  pro¬ 
viding  an  overall  summary  of  research  activities,  a  summary  of  key  findings,  and 
recommendations  for  subsequent  research. 

7.1  Research  Summary 

This  dissertation  provides  three  new  research  contributions.  First,  this  research 
demonstrated  an  appropriate  statistical  model  for  the  transformed  or  processed  pupil 
plane  data.  Pupil  plane  laser  speckle  is  transformed  by  computer  processing  to  pro¬ 
duce  images  related  to  noisy  object  autocorrelations.  The  statistics  of  the  non-imaged 
laser  speckle  is  well  studied  and  understood;  however,  the  Fourier-transform,  magni¬ 
tude  squared  operation  produces  the  noisy  autocorrelation  data  with  different  statis¬ 
tics.  This  document  provides  an  appropriate  mathematical  model  and  associated 
analysis  demonstrating  the  exponential  distribution  either  exactly  describes  or  well 
approximates  the  statistics  of  the  two-dimensional,  transformed  data  image.  This 
result  has  not  been  published  prior  to  this  research.  This  is  an  important  result  if  a 
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statistical  based  estimation  effort  is  to  be  considered  for  this  particular  phase  retrieval 
problem.  Also,  the  exponential  noise  model  has  implications  for  system  design,  pri¬ 
marily  related  to  SNR.  With  each  observed  data  frame  having  an  SNR  approximately 
equal  to  one,  SNR  is  improved  most  directly  by  increasing  observations  requiring  more 
operational  collection  time  required.  The  negative  exponential  noise  model  for  the 
noisy  autocorrelation  data  is  used  throughout  this  dissertation  to  include  producing 
data  for  simulation  efforts.  An  iterative  algorithm  using  the  exponential  model  was 
explored  without  success  in  Chapter  HU  However,  the  model  was  successfully  used 
throughout  the  research  to  include  a  theoretical  lower  bound  on  resolution  for  phase 
retrieval  from  correlography  data. 

Secondly,  this  research  provides  a  new  correlography  and  phase  retrieval  method 
using  polarization  diversity.  Image  recovery  using  polarization  diversity  has  been  ac¬ 
complished  previously  with  traditional  imaging  systems.  Employing  multi-channel 
diversity  in  pupil  plane  imaging  and  correlography  systems  has  been  previously  sug¬ 
gested  [23,24];  however,  this  research  demonstrated  a  novel  approach.  Also,  a  case 
was  made  here  for  a  simpler  detection  system  compared  with  previously  suggested  sys¬ 
tems.  Only  direct-detected  intensity  measurements  are  required  rather  than  difficult 
held  cross  products.  Using  active  illumination  with  polarized  laser-light,  man-made 
object  scenes  can  produce  polarization  diverse  reflections  and  scattering.  Observing 
these  effects  provides  additional  information  about  the  object  scene  enabling  improved 
estimation  efforts  in  the  ill-posed  phase  retrieval  problem.  Adding  a  secondary  ob¬ 
servation  channel  with  polarization  sensitivity  demonstrates  improvement  over  single 
channel,  polarization  insensitive  schemes.  Two  different  collection  schemes  were  in¬ 
vestigated  for  polarization  sensitivity. 

The  two  collection  schemes  motivated  two  iterative  algorithms  developed  us¬ 
ing  an  EM  approach.  A  new  technique  for  correlography,  borrowed  from  traditional 
image  recovery  [43] ,  used  a  polarization  parameter  to  relate  polarimetric  data  to  non¬ 
polarized  data.  This  added,  unknown  parameter,  has  the  potential  for  improving  the 
ratio  of  equations  to  unknown  variables  because  the  relationship  between  the  object 
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and  the  polarized  object  data  is  known  via  the  chosen  mathematical  model.  The  EM 
algorithm  method  provides  a  maximum-likelihood  solution  with  the  associated  con¬ 
vergence  properties.  Two  separate  EM  algorithms  were  presented  and  investigated 
with  simulation  and  laboratory  data.  The  first  algorithm  coupled  unpolarized  channel 
data  and  polarized  channel  data.  The  first  algorithm  employed  a  prior  distribution 
for  the  polarization  ratio  to  aid  the  solution  for  the  unknown  object  estimate.  The 
second  algorithm  employed  two  channels  both  sensitive  to  polarization  but  related 
by  90°  or  orthogonal  in  orientation.  A  previously  published  idea  for  algorithm  stop¬ 
ping  [31]  was  demonstrated  to  successfully  stop  both  algorithms  near  an  optimum 
iteration  number.  A  laboratory  experiment  was  used  to  demonstrate  the  validity  of 
approach.  Experimental  results  were  obtained  from  a  basic  laboratory  setup  for  the 
first,  two-channel  algorithm.  The  phase-retrieval  problem  still  remains  ill-posed  and 
computationally  difficult;  however,  adding  polarization  sensitivity  can  provide  sig¬ 
nificant  (greater  than  10%)  improvement  over  polarization  insensitive  schemes.  An 
operational  system  may  be  developed  with  available,  off-the-shelf  technology. 

The  third  research  contribution,  a  lower  bound  on  resolution  was  demonstrated 
for  phase  retrieval  with  both  polarization  sensitive  and  insensitive  systems.  Previous 
to  this  research,  bounds  on  recovered  image  resolution  have  been  demonstrated  for 
traditional  imaging  systems  obscured  by  atmospheric  distortion;  however,  this  has 
not  previously  been  done  for  non-imaging,  pupil  plane  systems.  A  lower  bound  on 
object  intensity  error  for  phase  retrieval  has  been  demonstrated  for  a  generic,  Gaus¬ 
sian  model  [5];  however,  this  research  used  a  more  accurate  statistical  model  and 
provides  resolution  as  a  measure  of  goodness  for  prospective  phase  retrieval  systems. 
The  computed  lower  bound  on  resolution  enables  comparison  between  polarization 
sensitive  and  insensitive  systems.  Comparing  the  computed  lower  bound  on  resolu¬ 
tion  for  a  simple  object  model,  this  research  demonstrated  theoretical  improvement 
with  the  introduction  of  polarization  sensitive  channel(s).  Statistical  independence 
or  statistical  diversity  is  of  primary  concern  for  multiple  data  channels. 
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7.2  Recommendations  for  Future  Research 

The  addition  of  polarization  diversity  is  shown  to  improve  the  phase  retrieval 
problem.  Additional  research  efforts  may  further  improve  this  difficult  problem  lead¬ 
ing  to  a  cost  effective  system  for  imaging  remote  satellites  at  geosynchronous  orbit. 
Suggested  follow-on  efforts  include: 

•  Investigate  additional  object  scene  diversity  such  as  range.  Range  diversity 
data  would  require  a  three-dimensional  LADAR  system.  Data  collection  in  a 
third  dimension  may  possibly  transform  this  effort  to  a  well-posed,  or  even  an 
over-determined  problem. 

•  Investigate  a  three-channel  polarimetric  system  in  hopes  of  improving  perfor¬ 
mance  with  additional  statistical  diversity.  A  three  channel  system  may  include 
an  unpolarized  data  channel  coupled  with  two  orthogonally  polarized  data  chan¬ 
nels. 

•  Further  investigation  into  solving  the  conditional  expectation  step  in  the  EM 
process  for  the  exponential  distribution  is  warranted.  Providing  an  EM  algo¬ 
rithm  with  the  exponential  density  function  may  provide  additional  improve¬ 
ment  over  the  Poisson  model. 
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Appendix  A.  Proof  of  Equation  \B. 111! 

The  development  found  in  this  appendix  is  primarily  from  [4],  This  appendix  pro¬ 
vides  the  proof  of  Eqn.  I2.ll  A  question  arises:  Does  averaging  the  lion-imaged 
laser  speckle  intensity  data  over  many  independent  realizations  achieve 
useful  information  about  the  object?  This  development  demonstrates  that  no 
useful  information  is  obtained  by  averaging  nonimaged  laser  speckle  patterns.  This 
development  ignores  photon  noise,  detector  read  noise,  and  dark  current  noise  dealing 
only  with  speckled  intensity  data. 

Theorem: 

K 

lim  K~l  //t(u)  —  C  (a  constant)  (A.l) 

X— >oo  ' 

k= 1 

where  K  is  the  number  of  collected  frames,  u  is  the  two-dimensional  coordinate  vector 
in  the  observation  or  pupil  plane,  and  p  is  kth  frame  of  the  observed  laser-speckle 
intensity  or  pupil-plane  data. 

Proof: 

As  previously  stated,  the  intensity  at  the  observation  plane  is  described  by 

/(n)  =  |  {/(x)}  |2  (A. 2) 

and  the  complex  field  at  the  object  plane,  /(x),  is  modelled  by 

/(x)  =  a(x)e^(x).  (A. 3) 

where  0  is  random  object  phase  uniformly  distributed,  U  ~  [ — 7r,  7t]  and  x  is  a  two- 
dimensional  coordinate  vector  in  the  observation  or  pupil  plane.  With  the  scaled 
Fourier  Transform  for  Fraunhofer  propagation  and  IA.31  the  field  at  the  pupil  plane  is 
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(A.4) 


n  u)  = 


a(x.)ej^]e-j2™x/Xz  rfx, 


where  A  is  the  optical  wavelength  of  the  illuminating  field  and  £  is  the  perpendicu¬ 
lar  propagation  distance.  Taking  the  magnitude  squared  of  Eqn.  |A.4|  produces  the 
observation  plane  intensity: 


J(u) 


a (Xl ) a (x2 ) ei0(xi } e“^(x2) e~j27ru(xi  ~X2 )/Xz  dxj dx2 , 


(A.5) 


where  xx  and  x2  are  independently  indexed  two-dimensional  coordinate  vectors.  To 
satisfy  Equation  1A,1|  the  ensemble  average  of  Equation  A,5i  is  taken  and  with  a  large 
number  of  independent  laser  speckle  patterns  is  equal  to  the  Expected  Value  operator. 


K 

lim  K~l  V  4(u)  =  E[  4(u)  ]  (A. 6) 

K — KX)  L ' 

k= 1 

It  is  assumed  independent  speckle  realizations  are  produced  by  minute  changes  in  the 
angle  or  phase  of  the  illumination  for  each  laser  pulse.  Taking  the  expected  value  of 
the  intensity  yields  the  following: 


E[I(  u)] 


E 


a(x1)a(x2)e^(xi)e-^(x2)e-j2,ru(xi“X2)/AVx1dx2 


(A.7) 


Using  linearity,  the  expected  value  operator  is  carried  inside  the  double  integral  oper¬ 
ation.  Also,  the  amplitude  of  the  reflected  field,  a  is  not  random  and  is  not  included 
within  the  Expectation  operator. 


POO  /»oo 


B[/(u)]  = 


a(Xl)a(x2) E  [eJ*(»>e-*,<®>]  e-J2™(*.-w)A*,jxi(fx2  (A.8) 


'  — oo  J  — oo 
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Analyzing  just  the  expected  value  operation,  this  is  easily  computed.  It  is  assumed 
the  surface  roughness  at  one  point  on  the  object  surface  is  independent  from  every 
other  point  on  the  object  surface.  Because  of  this  assumption,  </>(xi)  and  </>(x2)  are 
statistically  independent  when  Xi  ^  X2.  The  Expected  Value  Operation  becomes: 

£[e^(xi)]  •  E[e-^(X2)],  when  Xl  ^  x2  ( Case  1)  (A.9) 


£[eiW*,tW(x,f))]  =  E{ej0]  =  when  Xl  =  x2  ( Case  2).  (A.  10) 

Further  analysis  on  Case  1  above  and  the  use  of  Euler’s  identity  yields  the  following 
relationship: 


E[e^]  =  E[cos(<f>)  +  j  sin(0)]  =  E[cos((j))]  +  jE[sm((f))].  (A. 11) 


Finally,  assuming  the  random  surface  roughness  of  the  object  imparts  a  uniformly 
distributed  phase,  </>(x)  ~  U(— 7r,  7r) ,  we  can  compute  the  result  of  Case  1. 


E'[cos((/))] 


cos(0)  d(j) 


0 


(A. 12) 


i?[sin(0)] 


sin(0)  d(f> 


0 


(A. 13) 


Therefore,  Case  1  results  in  zero  and  Case  2  results  in  exactly  one.  This  can  be 
denoted  with  the  Dirac  delta  function,  5.  Substituting  the  Expected  Value  result 
back  into  Equation  IA.8  results  in 
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/oo  poo 

/  5(Xl  -  X2)a(x1)a(x2)e-^U(X1-X2)/A2rfx1dx2.  (A.14) 

oo  J  — OO 


Using  the  sifting  property  of  the  delta  function  simplifies  Equation  |A.  141  to  the  fol¬ 
lowing: 


E[I{  u)] 


a2(x1)e_j27ru  ('0')  dXl 
a2(x!)  dxx. 


(A.15) 


Analyzing  the  Equation  A. 151  the  integral  is  carried  out  only  over  the  hnite  region  i, 
illuminated  by  the  laser  beam  (this  is  the  beam  limited  scenario): 


a2(x)  dx  =  C.  (A.  16) 

C  is  a  constant  equal  to  the  object  intensity  or  brightness  function  summed  up  within 
the  region  of  laser  illumination.  Therefore  concluding  the  proof,  it  is  determined  that 
Eqn.  LA.  11  is  true. 


K 

lim  Ji-iy^/fc(u)=  C ,  (a  constant)  (A. 17) 

K— >oc  -i- — J 

k=  1 

From  this  proof  of  Equation  12.11  it  is  determined  that  averaging  nonimaged  laser 
speckle  data  collected  at  the  pupil  plane  over  many  pulses  yields  a  constant  value. 
The  processed  data  will  report  no  information  related  to  the  remote  object  as  the 
resulting  image  will  be  completely  washed  out. 

With  an  ergodicity  assumption,  this  result  is  also  achieved  if  one  performs 
a  time  average  of  a  single  speckled  intensity  collected  over  a  very  long  integration 


113 


period.  This  would  be  achieved  with  a  long  exposure  detection  and  illumination  with 
a  Continuous  Wave  laser  with  a  very  short  coherence  time. 
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Appendix  B.  Proof  of  Equation  \2.2 

This  appendix  provides  the  proof  of  Eqn.  12.21  This  result  first  appeared  in  [25] 
with  a  formal  proof  documented  in  [31].  Schulz  also  briefly  develops  this  result  [39]. 
Both  [31]  and  [39]  are  missing  a  term  in  the  solution;  therefore,  the  constants  are 
incorrectly  estimated.  The  published  results  also  ignore  the  photon  noise  generated 
during  detection  and  is  included  in  this  appendix.  The  development  found  here  is 
primarily  from  Phillips  [31].  This  proof  includes  photon  noise,  provides  the  missing 
mathematical  steps  not  found  in  Ref.  [31]  and  includes  the  missing  term. 

Theorem: 


K 

lim  l^_1-Uofc(u)}|2  ]  =  b\h(x)\2  +  c[R0(x)}  *  |/i(x)|2  (B.l) 

K — KX)  * 

k=  1 

JP-1  is  the  inverse  Digital  Fourier  Transform  (DFT)  that  is  performed  in  the  computer 
on  each  realization  of  the  observed,  speckled  Intensity  data,  IQ.  With  a  large  number 
of  laser  speckle  patterns,  we  replace  the  average  operation  with  the  expected  value 
operator. 


K 

lim  I<~‘  £[  |Jc~1{/„(u)}|2  ]  =  E\  |Jc“‘{/„(u)}|2  ]  (B.2) 

k= 1 

Proof: 

The  electric  field  at  the  object  plane  is  described  as 

/(x,t)  =  a(x,t)  •  e^(x)  (B.3) 

where  t  is  the  time  parameter,  x  =  (x,y)  is  a  two  dimensional  spatial  vector  in 
the  object  plane,  </>(x)  is  the  phase  directly  related  to  the  object  surface,  and  a(x) 
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is  the  field  amplitude.  We  are  assuming  the  field  amplitude  is  unknown  but  not 
random;  however,  the  phase  is  a  random  function  dependent  upon  the  object  surface 
height  profile.  The  phase  produced  by  each  laser  pulse  is  assumed  to  be  statistically 
independent  pulse  to  pulse  due  to  minute  changes  in  the  laser  illuminating  source 
(angle,  position,  etc.)  and  environment  (atmosphere). 

The  instantaneous  intensity  at  the  observation  plane  without  noise  is  described 
by 


I( u)  =  I  ^{/(x)}  |2  (B.4) 

where  TXz  is  a  continuous,  scaled  Fourier  Transform  representing  Fraunhofer  propa¬ 
gation  of  the  object  field  to  the  observation  plane.  The  observed  intensity  (measured 
through  optical  detection  devices)  is  formed  by  a  doubly  stochastic  process  distributed 
as  negative  binomial  (see  Section  [33).  The  observed  Intensity  data  is  modeled  as 

J0(u)  =  [/( u)  +  n(u)]  •  3(u)  (B.5) 

where  3(u)  is  the  aperture  function  denoting  the  region  where  the  speckle  pattern  is 
physically  recorded;  3(u)  =  1  for  the  points  within  the  measurement  aperture  and 
3(u)  =  0  elsewhere.  This  is  a  real  function  with  no  complex  phase  term.  Also,  n(u) 
represents  photon  or  shot  noise,  a  zero  mean  noise  such  that  the  observed  intensity 
(conditioned  on  the  average  photon  values)  has  a  Poisson  distribution  with  a  mean 
equal  to  the  intensity  without  photon  noise.  The  probability  density  function  of  the 
noise  function,  n,  is  unknown.  However,  the  noise  is  caused  by  the  random  arrival  of 
the  photons  as  they  emerge  from  the  optical  wavefront  onto  the  detector  device.  The 
previously  published  work  on  proving  Eqn.  |2.2  does  not  include  this  noise  term.  We 
include  it  here  for  completeness.  This  development  does  ignore  other  system  noise 
such  as  pre-amplifier  (read)  noise  and  dark  current  noise. 
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E{  |^-1{J0(u)}|2  ]  =  E[\  ^_1{-4(u)[/(u) +n(u)]}  |2] 
=  E[  Jr_i{A(ui)[/(ui)  +n(ui)]} 
x  «F{A(u2)[/*(u2)  +  n*(u2)]}  ] 


(B.6) 


Taking  the  DFT,  multiplying  the  two  terms,  and  distributing  the  expected  value 
Operator  leads  to  the  next  equation. 


N- 1  JV-1 

E[  |^-1{/o(u)}|2  ]  =  X  X  A(ui)A(u2)£?[/(ui)/*(u2)]exp{j27rw(ui 

ui=0 U2=0 
N- 1  N-l 

+  X  X  ^(ui)^4(u2)E[/(ui)n*(u2)]  exp{j27rw(ui 

ui=0  U2=0 
N-l  N-l 


+  X  X  ^(ui)^(u2)£[^*(u2Mui)]  exp{j27rw(u1 

ui=0  U2=0 
N-l  N-l 


+  X  X  ^(Ul)^(U2)^Kul)™*(u2)]  exp{j27rw(ui 

ui=0  U2=0 


u2)/N} 
u2)/N} 
u  2)/N} 

U  2)/N}  (B.7) 


We  will  hrst  investigate  the  last  three  terms  of  Eqn.  IB. 71  involving  noise.  The  cross 
terms  are  equal  and  exactly  zero  because  the  noise,  n(u)  is  independent  of  the  inten¬ 
sity,  J(u)  and  the  noise  is  zero  mean.  The  noise,  n(u),  is  generated  by  the  random 
arrival  of  the  photons  at  the  detector  surface  and  the  intensity,  /( u)  is  random  due 
to  the  uniformly  random  height  profile  of  the  object  surface.  Because  of  this,  we  can 
state  unequivocally  the  photon  noise  is  independent  of  the  random  phase. 


£[/(UlK(u2)]  =  E[I(u1)]E[n*(u2)]=0  (B.8) 

£[r(u2)n(Ul)]  =  E[I*{u2)]E[n(  u,)]  =  0  (B.9) 
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The  last  term  of  Eqn.  IB. 71  involves  the  autocorrelation  of  the  noise.  Assuming  the 
noise  at  each  detector  element  is  statistically  independent  from  the  noise  at  all  other 
detector  elements,  leads  to  the  result  in  Eqn.  IB.  101  where.  6  is  the  Dirac  delta  function 
and  a2n  is  the  noise  power  or  noise  variance. 


£[n(ui)n*(u2)]  =  cr2n  •  <5(ui)  (B.10) 


Substituting  the  result  in  Eqn.  IB.  101  back  into  the  last  term  of  Eqn.  IB. 71  produces  the 
contribution  due  to  photon  noise.  The  sifting  property  of  the  Dirac  delta  function 
simplifies  the  equation.  Also,  the  Fourier  Transform  of  the  aperture  function  is  a 
Dirac  function  at  the  origin. 


TV- 1  TV- 1 

EE  A(ui)A(u2)cr^(ui)  exp(j27rw(ui  -  u2)/N) 

ui=0  U2=0 

TV— 1 

=  an  ^(U2)  exp(-j27TWU2/A0 
u2=0 

=  ^(U) 


(B.ll) 


The  above  simplifications  of  the  noise  terms  simplifies  Eqn.  IB. 71  to  the  following 
equation. 


TV-1  TV-1 

E[  |.E~1{/0(u)}|2  ]  =  5^  5Z  ^(ui)^(u2)^[^(ui)/*(u2)]exp(j27rw(u1  -  u2)/iV) 

ui=0  U2=0 

+  <#(»)  (B.12) 


Most  of  the  previously  published  work  ignores  the  added  noise  term.  To  com¬ 
plete  our  analysis  of  the  Idell  function,  the  noise  term  is  temporarily  put  aside  to 


118 


compute  the  first  term  without  noise.  We  will  first  compute  the  autocorrelation  of 
the  speckle  intensity  without  photon  noise. 


Ri(u)  =  £[/(Ul)/*(  u2)]  (B.13) 

We  will  begin  expanding  the  equation  by  analyzing  the  intensity  at  the  observation 
plane. 


/(u)  = 


a(x)  •  e^(x) 


o-3 2^1 


dx 


(B.14) 


where  u  =  (w,  v )  is  a  two  dimensional  spatial  vector  in  the  observation  plane,  A 
is  the  wavelength,  and  2  is  the  propagation  distance  from  object  to  the  observation 
plane.  Expanding  Equation  IB.  131  enables  one  to  exchange  the  order  of  integration  and 
distribute  the  Expected  value  operator  inside  the  Fourier  integrals.  Only  the  phase 
term  of  the  reflected  field  is  random  and  subject  to  the  expected  value  operator. 


Ri(u ) 


xe 


-j2ir 


e^Tz1  P/2n 


clu2 


-3^ 


dxidx2dx1(ix2  (B.15) 


Analyzing  just  the  expected  value  operation  within  Equation  IB. 151  yields  the  term  T 
defined  as: 


T(Xl, x2, xi, x'2)  =  E[  e^(xi)e-^(x2 ]_  (B.16) 

There  are  5  cases  for  evaluating  T. 

1.  Xi  ^  x2  ^  x'  ^  x2 

2.  Xi  =  x2  =  x^  =  x2 
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3. 

Xl 

=  x2  and  x2 

=  x1  and  Xi 

x2 

4. 

Xl 

=  x)  and  x2 

=  x2  and  Xi 

x2 

5. 

Xl 

=  x2  and  x) 

=  x2  and  Xi 

~h 

/ 

xi 

Case  1.  If  all  spatial  vectors  are  not  equal,  then  we  assume  the  phase  terms  are 
statistically  independent.  Because  the  four  phase  terms  are  independent,  'h  becomes 
Eqn.  IB. 171  Because  the  phase,  0,  of  this  random  process  is  assumed  to  be  uniformly 
distributed,  ~  E[eJC<^]  is  exactly  zero  (c  is  any  integer). 


4/ !  =  E[  ei0(xi)  ]  ■  E[  e"^(x2)  ]  •  E[  }  ■  E[  ej^  ]  =  0  (B.17) 

Case  2.  If  all  spatial  vectors  are  equal,  then  the  exponents  add  to  zero  and 
the  expected  value  operator  yields  one,  denoted  with  the  Dirac  delta  function  for  the 
specific  condition. 


^  e(j?i(x1)-i0(xi)-j0(xi)+^(xi))  j 

=  E[ej'°] 

=  5(xi  =  x2  =  xl  =  x2)  (B.18) 

Case  3.  For  the  spatial  vectors  that  are  equal,  the  exponents  add  together. 
For  the  spatial  vectors  not  equal,  the  expected  value  operator  can  be  factored  into 
a  product  of  two  separate  expected  value  operations  due  to  statistical  independence. 
Again,  E[ec^]  =  0. 
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Vj/3  _  eO'2<^(x1)-j20(x2))  j 

=  £[ei2<Wxi)]  •  E[e~j2(t>ix2)} 

=  0  (B.19) 

Case  4.  This  case  is  similar  to  Case  3  as  it  factors  into  two  separate  expected 
value  operations;  however,  the  result  of  the  expected  value  operator  is  one,  as  in  Case 
2.  This  is  also  denoted  with  the  Dirac  delta  function;  however,  we  must  be  careful  to 
ensure  we  do  not  mathematically  duplicate  Case  2. 


\]>4  _  ^  e(j0(xi)-j</i(x2)-j</)(xi)+jXx2))  j 

=  E[ej-°]  ■  E[ej-°] 

=  S(xi  —  x[,  x2  —  x2),  (xiy^x2)  (B.20) 

Case  5.  This  case  is  the  same  as  Case  4.  This  is  also  denoted  with  the  Dirac 
delta  function.  We  must  be  careful  to  ensure  no  duplication  with  Case  2. 


^  ^  eO'</>(xi)-j</>(xi)-jy(x'1)+j0(x,1))  j 

=  E[ej  0]  ■  E[ej'°] 

=  S(x !  -  x2,  x[  -  x'2),  (x1^x'1)  (B.21) 

Non-zero  results  are  obtained  from  Cases  2,  4,  and  5.  Cases  1  and  2  do  not 
contribute  to  the  equation.  Substituting  the  Dirac  delta  functions  back  into  Equation 
B.15lwill  enable  some  simplifications.  The  sifting  property  of  the  Dirac  delta  function 
will  reduce  variables  and  integrands. 
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poo  poo 


Ri{  u)  = 


a2(x!)a2(x2)e 


-J'27rxi  (UiazU2)  eJ27rx2  (U\,U2) 


dXiC?X2 


r  — OO  — OO 

poo  /»oo 


+ 


a2(xi)a2(x1)rfxirfx1  — 


a4(xi)dxi 


(B.22) 


'  —  OO  «/  — OO 


The  first  two  terms  in  Equation  IB. 221  are  from  Cases  4  and  5,  respectively.  The  third 
term  is  from  Case  2  and  must  be  subtracted  once  because  it  appears  as  a  special 
case  in  both  the  first  and  second  terms.  The  results  from  Cases  2,  4,  and  5  are  now 
accounted  for  without  duplication.  Note,  Term  3  is  the  missing  term  not  found  in 
previously  published  papers.  Now  that  we  have  Equation  B.22,  the  autocorrelation  of 
the  intensity,  we  substitute  this  result  back  into  Eqn.  IB. 121  We  will  ignore  the  noise 
term  temporarily  and  add  it  back  in  at  the  end  of  this  proof.  Each  term  is  dealt  with 
separately. 


N—l  N— 1  ^oo  „qq 

e\  |J-'{/„(u)}|2  ]='£'£  / 

,, _ n  ,, _ n  J  —  oo  J  —  oo 


a2(xi)a2(x2) 
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Term  1.  We  will  simplify  the  integral  with  substitution  of  variables.  If  we  let 
Xq  =  X\  —  x2,  Term  1  becomes 
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where  i?o(x0)  is  considered  the  autocorrelation  of  the  object  intensity  with  a  spatial 
difference,  xq.  We  will  continue  to  group  terms  and  simplify. 
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(B.25) 


Analyzing  just  the  term  in  brackets,  it  is  a  scaled  Fourier  Transform  of  the  intensity 
autocorrelation  which  is  defined  as  |F0(£)|J.  Note,  since  Ra  is  a  symmetric  function 
about  the  origin,  its  Fourier  Transform  is  symmetric. 
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(B.26) 


Plugging  this  result  back  into  Term  1  yields 
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Using  variable  substitution,  let  u0  =  ui  —  u2. 


N- 1  N- 1 
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ui=0  U2=0 

Next,  group  terms  as  we  have  previously  done. 
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where  iUi(uo)  is  the  autocorrelation  of  the  aperture  function  with  a  spatial  difference 
of  Uo-  If  we  ignore  the  normalization  constant,  this  is  recognized  as  the  diffraction- 
limited  Optical  Transfer  Function  (OTF),  7i.  Also,  after  ignoring  the  normalization 
constant,  the  inverse  Fourier  Transform  of  the  OTF  is  the  Point  Spread  Function 
(PSF)  or  impulse  response  multiplied  by  a  constant  [18]. 
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|/i(w)|J  =  7- — --  /  Ra(uq)  exp(— j— wu0)du0  (B.31) 

(A z)A  J \z 


R-'iHiv)}  =  (A z)2  \h(Xzw)\2  (B.32) 

Finally,  for  Term  1,  we  recognize  that  Equation  IB. 291  is  an  inverse  Fourier  Trans¬ 
form  of  the  product  of  two  symmetric  functions.  Applying  symmetry  and  the  convolu¬ 
tion  theorem  of  the  Fourier  Transform  yields  Eqn.  B.331  We  have  the  autocorrelation 
of  the  object  brightness  function  imbedded  within  our  result. 
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T\  =  (Xz)2R0  (Acw)  *  \h(Xzw)\2 

=  ci?0(x)  *  |/i(x)|2  (B.33) 

Here  in  this  result,  c  =  (A^)2  and  x  =  Xzw.  Therefore,  Term  1  is  a  constant  factor 
times  the  impulse  response  convolved  with  the  scaled  autocorrelation  of  the  object 
intensity. 

Term  2.  We  can  simplify  the  integration  by  using  another  substitution  of 
variables,  similar  to  what  was  done  for  Term  1.  Let  xq  =  xi  —  x[. 
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Again,  Ro(x0)  is  the  autocorrelation  of  the  object  intensity  with  spatial  difference, 
xq.  Using  another  substitution  of  variables,  let  uq  =  Ui  —  u2  and  group  terms. 
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where  -Ra(uo)  is  the  autocorrelation  of  the  aperture  function  with  a  spatial  difference 
of  u0.  Again,  this  is  recognized  as  the  diffraction-limited  Optical  Transfer  Function 
(OTF),  7i.  ignoring  the  normalization  factor.  As  stated  previously,  the  inverse  Fourier 
transform  of  the  OTF  is  the  scaled  impulse  response.  Next  we  will  define  Sr  as  the 
total  sum  of  the  object  intensity  autocorrelation.  This  is  a  constant  value  since  the 
object  intensity  is  finite  in  extent  due  to  the  beam  limited  scenario. 


Sr  = 


-Ro(x0)dx0 


(B.36) 


Finally,  for  Term  2,  we  are  able  to  simplify  to  an  easily  described  equation. 


T2  =  SR-  (\z)2\h(\zw)\~ 
=  SR  •  c\h  (x)  j 2 


(B.37) 


Term  3.  Term  3  is  the  missing  term  not  found  in  published  literature.  We  will 
simplify  Term  3  by  first  defining  a  term  of  summed  object  intensity  squared,  Sp. 


Sp 


(B.38) 


This  is  a  constant  term  and  can  be  factored  outside  any  remaining  integrands.  We 
realize  this  quantity  is  finite  since  we  are  in  a  beam  limited  scenario.  Returning  to 
Term  3  and  rearranging  terms  to  identify  the  summed  intensity  squared  term  we  find 
the  following: 
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We  will  continue  the  simplification  with  a  substitution  of  variables  by  letting 
u0  =  tp  —  u2  and  rearranging  terms. 
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u0=0 

Again,  .Ra(uo)  is  the  autocorrelation  of  the  aperture  function  with  a  spatial  difference 
of  Uo  or  the  diffraction-limited  OTF  without  the  normalization  factor.  Recognizing 
the  last  equation  is  an  inverse  Fourier  Transform  of  the  OTF,  Term  3  simplifies  to 


T3  =  —Sp  ■  ( Xz)2\h(Xzw)\ 2 
=  —  Sp  ■  c\h(x)\2. 


(B.41) 


The  constants  from  Terms  2  and  3  can  be  combined  into  a  single  constant,  b. 


b  =  c 


R0(x)dx  — 


(B.42) 
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Finalizing  the  Idell  Function.  Combining  the  three  terms  together  to  yield 
the  final  result  for  F'[|JF_1{J0(u)}|2]  (ignoring  photon  noise)  produces  the  final  equa¬ 
tion. 


E[ \T  1  (d0(u) } | ”]  =  b\h  (x) | 2  +  c[R0  (x)]  *  |h(x)|2  (B.43) 

Therefore,  the  published  result  is  proven.  However,  the  constant,  b ,  is  not  defined 
correctly  in  previous  literature  due  to  the  missing  term  defined  above. 

Adding  the  noise  term  due  to  photon  noise,  the  final  result  is  given  by 


E[ \T  1{I0{u)}\2]=b\h(x)\2  +  c  [R0  (x)]  *  |/i(x)|2  +  a25(x).  (B.44) 
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Appendix  C.  Proof  of  Exponential  Statistics  for  the  Processed  Data 

This  appendix  details  the  proof  for  the  transformed  laser  speckle  data  or  noisy  auto¬ 
correlations  distributed  as  exponential  with  a  mean  equal  to  the  autocorrelation  of  the 
true  object.  This  proof  assumes  each  laser  speckle  image  is  statistically  independent 
spatially. 

A  single  image  of  laser  speckle  data  is  approximated  very  well  to  be  distributed 
as  negative  binomial  [17].  In  most  detection  schemes,  the  speckle  will  exhibit  small, 
localized  correlation;  however,  for  this  proof  we  will  assume  that  each  point  in  the 
laser  speckle  image  is  statistically  independent.  To  form  the  autocorrelation  data, 
each  of  the  observed  laser  speckle  intensity  images  are  post-processed  by  DFT  and 
magnitude  squared  operations.  To  prove  the  statistics  of  the  processed  data,  this 
section  follows  Goodman’s  treatment  of  random  phasor  sums  [17]. 

Each  point  in  the  intensity  data  is  a  real-valued  but  random  number  (phase 
equal  to  zero).  The  Fourier  Transform  kernel  provides  a  complex  but  known  phasor. 
For  this  analysis  we  will  assume  the  two-dimensional  data  image  is  JV  x  A  with  N 
sufficiently  large  to  utilize  the  Central  Limit  Theorem.  The  sum  of  a  large  num¬ 
ber  of  independent  random  variables  (in  this  case,  the  random  intensity  values)  is 
asymptotically  Gaussian  as  N  — >  oo  [17]. 
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Equation  IC71  is  essentially  equivalent  to  a  random  phasor  sum.  The  intensity  IQ  is 
a  collection  of  random  variables  where  the  phase  is  known  (exactly  zero)  and  the 
magnitude  is  random  but  statistically  independent  of  each  other.  To  discover  the 
statistics  of  the  transformation  produced  by  Eqn.  IC.ll.  we  will  execute  Goodman’s 
random  phasor  sum  procedures  [17].  The  real  and  imaginary  parts  of  the  Fourier 
transformed  data  are  defined  as  r  and  i,  respectively. 
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N-l 

r(w)  =  Re[  E~l{I0}  ]  =  ^  I0(u )  cos 

u= 0 

N-l 

i(w )  =  /m[  jr_1{/0}  ]  =  ^  J0(w)  sin 

ii=0 

Next,  the  mean  of  the  real  and  imaginary  parts  are  computed.  Because  the 
expected  value  of  laser  speckle  is  a  constant  (e.g.  E[I0]  =  C ),  the  expected  value 
of  the  real  and  imaginary  parts,  are  two-dimensional  Cosine  and  Sine  transforms, 
respectively,  of  a  constant  value.  The  mean  of  the  real  and  imaginary  parts  are  equal 
except  at  w  —  0, 


27 TUW 
N 

27 TUW 
N 


(C.2) 

(C.3) 


E[r(w)}  =  CN2S(w),  (C.4) 

E[i(w)\  =  0  V  w.  (C.5) 

The  variance  of  the  real  and  imaginary  parts  are  similarly  computed.  The  au¬ 
tocorrelation  of  the  laser  speckle  is  approximated  by  a  delta  function  plus  a  constant, 
due  to  assumed  statistical  independence  and  the  second  moment  of  the  speckle  also 
equal  to  a  constant  (e.g.  E[Iq\  =  B ), 

E[I0(u)I0(v)\  =  C2  +  B5{u-v).  (C.6) 

The  variance  of  the  real  and  imaginary  parts  are  computed  to  be: 
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Therefore,  the  variances  of  the  real  and  imaginary  parts  are  equal  except  for  a  small 
number  of  points  in  the  image, 


var [r(w)]  =  va,r[i(w )]  Vw9  {0, iV/2}.  (C.9) 

Next,  the  correlation  coefficient,  pri,  for  the  real  and  imaginary  parts  is  com¬ 
puted  to  be  zero  everywhere  due  to  orthogonality  and  are  uncorrelated, 


TV- 1  TV-1 

pri(w )  =  E[r{u)i(y)\  =  EE  E  [I0(u)I0(v)\  cos 

It— 0  -u=0 

=  0.  (C.10) 

Therefore,  a  majority  of  the  Fourier  transformed  image  points  (V  w  3  {0,  iV/2}) 
are  jointly  distributed  as  circular  complex  Gaussian  random  variables.  After  the 
DFT  is  taken  on  the  laser  speckle  data,  a  second  operation,  magnitude  squared,  is 
performed  to  finalize  the  transformation  from  laser  speckle  to  noisy  autocorrelation 
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data.  The  magnitude  squared  operation  transforms  the  circular  complex  Gaussian 
random  variables  to  an  intensity  distribution.  The  intensity  distribution  of  circu¬ 
lar  complex  Gaussian  random  variables  obey  exponential  statistics  [17].  The  mean 
of  these  exponential  random  variables  is  the  autocorrelation  of  the  true  object  as 
discussed  in  Section  13.21  The  few  pixels  in  the  image  where  the  imaginary  part  dis¬ 
appears  is  less  than  ten  percent  of  the  image.  With  the  assumptions  made  above  and 
ignoring  the  center  pixel  where  the  peak  of  the  autocorrelation  occurs,  a  single  frame 
of  noisy  autocorrelation  data  can  be  well  approximated  as  exponentially  distributed 
random  variables. 
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Appendix  D.  Single  Channel,  EM  Phase  Retrieval  Algorithm 

This  appendix  details  the  single  channel,  EM  phase  retrieval  algorithm  using  Poisson 
statistics.  The  result  is  nearly  identical  to  the  ML  algorithm  presented  by  Schulz 
and  Snyder  [38]  only  differing  by  a  scale  factor.  This  algorithm  is  re-derived  here 
for  completeness  and  comparison  to  the  multi-channel  EM  algorithms  developed  in 
Chapter  [TV] 

The  complete  data,  dk(y,x),  is  chosen  such  that  it  is  related  to  the  incomplete 
data,  dk(x)  by 


dk(x)  =  ^dk(y,x),  (D.l) 

y 

where  k  denotes  the  data  frame,  x  and  y  are  two-dimensional  spatial  variables,  and 
the  complete  data  is  assumed  to  be  independent  and  identically  distributed  Poisson 
random  variables.  Knowing  information  about  the  observed,  incomplete  data,  we 
choose  the  mean  of  the  complete  data  to  be 


E[dk(y,  x)]  =  o{y)o(y  +  x).  (D.2) 

Because  the  complete  data  is  chosen  to  be  a  Poisson  random  variable  with 
a  known  mean,  we  can  completely  describe  the  probability  mass  function  of  the 
complete  data. 
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The  observed  incomplete  data  is  described  by  its  mean 
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(D.4) 


E[dk(x)]  =  y^E[4(y,i)] 
y 

=  ^2o(y)o(y  +  x) 
y 

The  log-likelihood  function  of  the  complete  data,  Lcd,  is  found  by  taking  the 
natural  log  of  the  probability  mass  function  in  Eqn.  ID. 31 


Lcd(o)  =  x )  loS  [°(y)o(y  +  a:)]  -  o{y)o{y  +  x)  +  A.T.  !•  (D.5) 

Expectation  Step.  The  expectation  step  of  the  EM  algorithm  is 
defined  as  the  expectation  of  the  complete  data  log-likelihood  function  conditioned 
on  the  old  object  estimate  (from  the  previous  iteration),  o°ld(y),  and  the  incomplete 
data,  dk(x). 


Q  =  E[Lcd  \  o°ld,dk(x)] 

=  S  {  E°M  x)]  loS  [°(y)°(y  +  x)]  “  o(y)o(y  +  x)  +  A.T.  j  (D.6) 

where  Eold  is  the  conditional  expectation  conditioned  on  the  old  estimates  of  o  and 
the  incomplete  data.  K  is  the  total  number  of  frames  and  A.T.  denotes  another  term 
not  a  function  of  o. 

For  Poisson  statistics,  the  conditional  expectation  on  the  complete  data  given 
the  incomplete  data  is  detailed  by  Shepp  and  Vardi  [42]: 
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Maximization  Step.  The  maximization  step  involves  maximizing 
the  Q  function  with  respect  to  the  unknown  parameter,  o.  We  take  the  partial 
derivative  of  Q(o|o°ld),  set  it  equal  to  zero  and  solve  for  the  object,  o. 
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where  S “ew  is  the  sum  of  the  new  object  estimate  defined  as 
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Setting  Eqn.  ID. 81  equal  to  zero  and  solving  for  the  object,  o,  yields 
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In  order  to  simplify  this  expression,  a  few  new  terms  will  be  defined.  We  will 
let  the  average  of  the  incomplete  data  be  defined  as  the  measured  autocorrelation, 
Ro, 


R°(x)  =  jr  y^4(V)-  (D.ll) 

k=  1 

The  autocorrelation  of  the  old  object  estimate  from  the  previous  iteration  will  be 
defined  as 


RT(x)  =  J2o0ld(y)o0'd(y+x)-  (D-12) 

y 

Also,  we  will  define  the  following  parameter 
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Substituting  these  parameters  back  into  the  expression  for  the  object  estimate 
produces  the  following  equation. 
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Here,  the  expression  for  the  new  object  estimate,  onew,  is  a  function  of  the  old 
object  estimate,  o°ld,  and  the  incomplete  data.  Summing  both  sides  of  Eqn.  ID, 141 
provides  an  expression  for 


Q  new  _ 


(D.15) 


The  iterative  solution  restated  is 


'(ifo)  = 


2S" 


0°*  *  2k  +  Oou  *  R' 


R° 


R°ld 


(yo 


(D.16) 


This  result,  is  identical  to  the  original  Schulz  and  Snyder  ML  algorithm  [38]  differing 
only  by  a  scaling  parameter,  S “ew.  The  Schulz  and  Snyder  ML  algorithm  uses  the  S'°ld 
scale  parameter.  The  Schulz  and  Snyder  ML  algorithm  is  detailed  in  Eqn.  ID. 171 


'(yo)  = 


2S°ld 


/nold  -x- 
O  ik 


Rq 

R°h 


+  oold  * 


Rq 

R°h 


(: Vo ) 
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